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THERMAL-STRESS ANALYSIS FOR A WOOD COMPOSITE BLADE 


Introduction 

Since 1977, NASA Lewis Research Center has been pursuing the development 
of low cost rotor blade technology for large horizontal axis wind turbine. 
Laminated wood manufactured by bonding together 1/10" to 1/8" thick sheets, or 
plys, with epoxy is a particularly attractive candidate material for a rotor 
blade since the raw material has high specific strength, high specific stiff- 
ness and low cost. This led to the manufacturing of eight laminated wood 
blades in 1980. They were installed as Mod-OA wind turbine at three different 
locations; Kahuku Point, Hawaii, Culebra, Puerto Rico, and Block Island, Rhode 
Island. Mod-OA turbines have 125 1 diameter rotors and generate 200 Kilo-Watts 
each when in operation. Their performances were closely monitored. In 1981, a 
crack was found in Blade No. 1012 located in Block Island, which occurred in 
the leading edge extending along the entire blade. The cause of the crack was 
unknown, but thermal-stress induced by the sun vias suspected to be one of the 
main reasons. Therefore, the investigation has been pursuaded along this 
line. 

The analysis of thermal -stress of the blade induced by solar insolation 
consists of two phases. The first phase is to find the temperature distri- 
bution throughout the blade, a heat conduction problem. The second phase 
is to determine the thermal -stress distribution of the blade caused by the 
temperature distribution found in the first phase, a thermal -stress analysis 
problem. Since wood is an orthotropic material in the senses of heat 
conduction and stress-strain relationships, these characteristics must be 
included in the analyses of both phases. 



Phase X: Heat Conduction Analysis 


A. The Blade and the Assumed Environment 

Blade No. 10X2 is 700 inches long with maximum chord width of 62.4 inches, 
see Fig. 1 for its overall dimensions. Airfoil standard section NACA 230 [1], 
is used. The thickness to chord ratio of the blade is varied from 31.75 at 
station 150 to 7.5 at station 732. 

In the determination of temperature distribution throughout the blade, it 
is reasonable to assume that a two-dimensional heat conduction analysis of a 
typical blade section is adequate, because the blade is basically a slender 
body. In addition, the thermal conductivity of wood in longitudinal direction 
of the blade, which is parallel to the grain, is two to three times higher than 
those in transversal directions. As a result, the temperature gradient in 
longitudinal direction is far smaller than those in transversal directions. 

This again indicates that a two-dimensional analysis is adequate. 

The section at station 126 is arbitrarily selected as a typical section 
for analysis. The thickness to chord ratio at this station is 30.4, which is 
close to those of the sections in the neighborhood of station 300, the main 
body of the blade, see Fig. 2. In the nose portion, called 0-spar, 3" thick 
laminated rotary veneers of Douglas Fir is used. In the tail portion, 3/4" 
thick honeycomb paper core is installed to provide stiffness and light-weight. 
Pieces of sawn Douglas Fir are used for stringers. On the surface, the blade 
Is covered with a layer of 1/8“ thick Birch Plywood. 

Two parked positions of the blade are considered; one keeps the section In 
horizontal position, the other, vertical. Technically speaking, the difference 
between these two positions is that the heat flux inpui: zones are different as 



shown in Fig. 3 and 4. The heat flux input *one is defined as the portion of 
boundary surface of the blade exposed to the sun. The solar insolation is 
estimated at 363 BTU/y^-hr [23, The absorbtivity of the blade surface is 
estimated at 0,9, which is not only on the safe side for the purpose of stress 
analysis but also has a good reason: after a long period of exposure, the 

blade surface is usually dirty and full of bug remains, which often drastically 
increase the surface absorbtivity. Thus, the heat flux input, q r , for the heat 
flux input zone is 


q r - 363 x 0.9 = 326 BTU/y^Z-^j, . 

The ambient temperature, T , is assumed to be 90°F, typical temperature for 

co 

a hot summer morning in Block Island. Further, the air is assumed to be 
still, arid the convective heat transfer coefficient, h, for the blade surface 
is assumed to be 5 BTU/f^2-f 1 r.°F * 

In description of wood properties, three mutually perpendicular axes, 
longitudina' „ radial and tangential, need to be used. The longitudinal axis, 

L, is parallel to the grain; the radial axis, R, is normal to the growth rings; 
and the tangential axis, T, is perpendicular to the grain but tangent to the 
growth rings. The thermal conductivities of the wood in L, T and R directions 
[9] are listed in Table 1, note that a small difference of 10% is assumed 
between ky and k^ to account for possible physical differential between these 
two directions. Since rotary veneer is used in blade construction, the axes 
L, T, and R correspond to principal material axes c and n respectively. 

The material axes will be used in next section. 

The thermal conductivity of honeycomb paper core may be estimated by 
using the average value for paper and air, which is 0.013 BTU/f^f^ y . 



Table 1. Thermal Conductivity of Wood 


Wood 

k L 

k T 

k R 

Douglas Fir 

0.172 BTU/ hr . ft> . F 

0.068 BTU/fop.ft. *F 

0.075 BTUV.ft.’p 

Birch 

0.211 

1 

0.084 

0.093 


B. Governing Partial Differential Equation 

Because of the configurational complexity of the blade section and the 
orthotropic property of heat conduction of the wood composite, finite element 
method with variational approach is chosen as the method of solution. This 
method may also be called as a numerical method in solving partial differential 
equations. For an element of an orthotropic material, the controlling partial 
differential equations may be derived in the following outlines. 

Heat flux input - Heat flux output = Heat retention in body (1) 

Referring to Figure 5, let ‘f l be the heat flux and £ k‘, the thermal 
conductivity, their directions along axes y, z, 5 and n are indicated by 
subscripts. 'T 1 represents the temperature, ‘Cp 1 , specific heat and *p', 
the density of wood. Thus, Equation (1) may be expressed as 

fj dxdz + f z [ dxdy - Cf v t dxdz] - Cf z | dxdy] 

j y z z * y+Ay 4 z+Az 

3T 

= p dxdydz C n — (2) 

p 3t 

Applying Taylor's Series and simplifying, Equation (2) becomes. 
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(3) 


Since the principal axes of the material, 5, and n, are in general at an angle, 
3, with global axes y and z as shown in Fig. 5, the following transformation 
may be used. 





Using chain rule; 


8T 


~*3y 

9z 


OT 






OT 

— — 





— 


cos 3 

sin 

3 



35 


35 

35 


ay 






ay 

OT 


ay 

Oz 


3T 


-sin 3 

cos 

3 


OT 

On 


on 

on 


Oz 






32 


we have. 


cos 3 -sin 3 
sin 3 cos 3 


*c 0 


0 -k 


n 


cos 3 sin 3 

■sin 3 cos 3 



il 


oy 


OT 


Oz 


(4) 
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Substitute (4) into (3), we obtain: 

3 2 T 3 2 T 

tkr cos 2 g + k n sin 2 e) — - + (k r - k^) Sin2B ----- 
5 1 0y2 ^ 1 3y 9z 


3T 

+ (kg sin 2 P + k n cos 2 £) --- = p C p -- 

^ U tc 


3 2 T 

3z 2 


(5) 


Equation (5} is the governing partial differential Equation for the heat 
conduction analysis of the blade. 

The boundary conditions around the surface of the blade are: 


(a) prescribed heat flux, q r , in the heat flux input zone, namely, 

q r = 326 BTU/ ft 2 hr . (5a) 

(b) convective heat transfer loss, q c , for the entire boundary 
surface of the blade, 

q c = h(T s - TJ (5b) 

where T s is the surface temperature of the blade, a variable. 


C. Calculus of Variation 

For a given problem with a governing differential equation and boundary 
conditions, the task of the variational formulation is to find the unknown 
function *F* for the differential equation, which extremizes or makes 
stationary a functional 'I 1 subject to the same boundary conditions of the 
problem. 
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where. 


I 


-t 


F(x, T(x), 


T‘{x)) dx 


T(x), is a weak variation of T(x}, namely 


( 6 ) 


T( Cx) - T{x) + en(x) (7) 

e is an error term of a small magnitude , n(x) is an error function which 
vanishes at its boundaries. 

From the calculus of variation [3], when 'I 1 is extremized respect 
to ! T* and the error term e is made to approach zero, the Euler-Lagrange 
Equation is formed, namely 


9F d 9F 

9T dx ai' 


( 8 ) 


In other words, the functional 1 1 * is extremized or made stationary when 
the Euler-Langrange equation and its boundary conditions are satisfied. 
Therefore, the task of seeking solution, T{x), for the governing partial 
differential Equation (5) becomes a numerically simpler problem, that is, 
finding first the 'F' in such a way that its Euler-Lagrange equation is 
identical to the governing differential equation. Next, the Functional ‘I* 
Is formed. Then the finite element method is used to find an approximate 
temperature profile which extremizes the functional and satisfies the 
boundary conditions. The temperature profile so obtained is the solution 
to the heat conduction problem in question. 



For the governing differential equation, (5), the function 'F 1 may be 
formed as follows: 

1 3T 2 3T 2 

F - - C(kr cos 2 3 + k n sin 2 3} {--) + (kr sin 2 3 + k n cos 2 3) (--) 

2 ^ 1 3y 3z 

81 8T 3T 

+ (kr " k n ) sin23 (--) (--} + 2p C d -- T] (9) 

^ 11 3y 8z y at 

Thus, the functional which contains the Function 'F 1 as described 
in Equation {9) and the corresponding function for boundary conditions 
(5a) and (5b), is 


I = / FdV + / [q r T + x h (T “ T ) 2 ] ds 


(10) 


D. Finite Element Method - Variational Approach 

The cross section of the blade at station 126 is selected as the typical 
section for temperature profile determination. Triangular elements are used. 
The blade consists of 3 different kind of materials, Douglas Fir, Birch 
plywood and honeycomb paper core. It was divided into 75 regions, see 
Fig. 6. Each region was further divided into proper numbers of columns and 
rows. See Table 2. Then the triangular finite element mesh is generated 
automatically [43. The results of regions 12 and 54 are shown in Fig. 7 as a 
typical example. 


Table 2 


9 
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To describe the finite element method,, a typical element as shov/n in Fig. 8 
is chosen. An approximate temperature profile, T(y,z), can be constructed by 
using shape functions Nj, Nj and M k to linearly interpolate among the nodal 
temperatures Tj, Tj and T k , the values of which are to be determined. T(y,z) 
is described as: 

T(y,z) * N 1 * Tj + Nj Tj + N k T k (11) 


where 


Mi = — [a n - + b.,* y + c* z] 
2A 1 1 1 

“5 = ll Ca j + b 5 y + c 3 ^ 


(11a) 


h s ^ Cajc + b k y + c k ^ 

1 Yj Zj 

1 Y j z j 
1 Y k Zk 


- Y A 

~ Y k Zj, 

a j = Y k Zj 

2 k Yj, 

a k = YjZj - YjZ 


“ z k . 

b j = z k 

- h , 

b k = Zj 

" z j 

= Y k 

Y j » 

Cj = Yj 

“ Y k , 

c k = Y j 

“ Y f 



(lib) 


(He) 


Substitute T(y,z) into Equation (10), and partial differentiate l I’ 
with respect to unknown nodal temperatures. Then, set them equal to 
zero. The result is a set of equations: 


CC] + CK] [T] = CP] 

dt 
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CC3 »f [C e ], [K3 = H K e ], 

e“l e-1 


CP] = I CP e ] 

e=l ' 


2 1 1 

pCnA 

[C e ] " 121 , 

1 1 2 


[P e ] = (-q r + h T } * 




cos 2 3 + kp sin 2 3 


b i b i b i b j b 1 b k 

b jb-f bjbj bjb k + 

— b k b i b k b j b k b k_ 


k ? sin 2 3 + k n cos 2 3 
4A 


” c i c i 

CjCj 

c i c k 

(kr - k n ) sin23 

bf c i 

b i G j 

b i c k ” 

c j c i 

c i c j 

c j c k 

+ 

4A 

bjc-j 

b j c j 

b j c k 

_ c k c i 

c k c j 

c k c k_ 


_ b kci 

b k c j 

b k c k _ 


h Ljj 


2 10 

12 0 or + 
0 0 0 


b _Ljk 

"i 


ooo" 
0 2 1 
0 12 


“2 0 r 
0 0 0 
10 2 


h L ki 
or + — — 
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Note that n is the total number of elements of the blade. [Cl is the 
capacitance matrix, [Pi, the heat flux matrix and [Kl, the conductance matrix 
of the entire blade. [P e l and [K e ] contain 'or 1 in their expressions. The 
purpose is to accommodate the elements with one edge or two edges on the 
boundary of the blade. For example, if edge j-k is on the boundary, the term 
which contains is maintained and the other two are deleted. For a non- 
boundary element, L-j-j, Lj^ and L^* are assumed a value of zero. Thus, the 
last term of [K 6 ] is deleted and [P e ] is dropped off for that element. 

Equation (12) describes the heat conduction problem in trancient state. 
When it reaches the steady state, the rate of change of temperature respect to 
time vanishes, i.e., 

9{T> n 

“ ft 

at 


Equation (12} becomes. 


[K3 [T] - [Pi 


Equation (13) describes the heat conduction problem in steady state. 


(13) 


E. Heat Conduction in Steady State 

Two steay-state Heat Conduction Analyses are made on the University NAS 
6650 computer system. The first analysis is performed when the blade section 
is placed in horizontal position, the second, in vertical position. The sun 
is assumed to shine vertically on the top of the blade, which defines the heat 
flux input zone as shown in Fig, 3 or 4. So only the element with an edge 
on the top side has heat flux input 1 q r * entering the element thru that edge. 
The temperatures on all nodal points were obtained. Only the temperatures of 
the exterior nodes are plotted in Fig. 9 and 10. Temperatures of interior 
nodes are shown representively on typical sections in Fig. 11. 
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Observing Fig. 9 and 10, one finds that both of the sunny side surfaces 
reach a high temperature of 155 “F and the opposite sides are on or slightly 
above 90“F. The identicalness of the results and that the low temperature is 
so close to the ambient temperature offer strong assurances that the analysis 
is correct. Further, Fig. 11 shows a smooth transition from high temperature 
side to low temperature side in all sections. They again indicate that the 
results are quite reasonable and hence reliable. 

F. Heat Conduction in Transient State 

It is important to determine the length of sunshine time required for the 
blade to reach the steady-state temperature. If the length of time required 
is not more than the sunshine hours of a day in the summer, the steady-state 
temperature may be used as the critical temperature distribution for anlaysis. 
Otherwise, temperature distribution at the end of sunshine hours of a day 
should be used. 

In time domain, a linear interpolation model is used, see Fig. 12. The 
unknown temperature field (T>, spanning from one time to another separated 
by a time step A t, may be expressed as 

(Tl = N 1 (T-j } + Nj lTj> (14) 

where N-j and Nj are shape functions, namely, 

t 3N* -i 

N,- = 1 - - --- = - 

1 At 3t At 

hence, 

fj = i. 

at At 
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substitute Equation (14) into Equation (12) and apply Galerkin's method [5], 
two integral matrix equations are produced, they are. 


At d(T} 

/ Ni ([C] ---- + DCHT) - {P» dt “ 0 
o dt 


(15a) 


At d{T} 

/ N n * ([C] f DGlT} - {P}) dt = 0 

o dt 


(15b) 


The time derivative of Equation (14) is. 


dtn dNi 

dt dt 


dM -f X 

--- TT-? > + --- rr,-} = — {Ti } + -- tT-f} 
dt J at 1 a+ J 


1 

At 


(16) 


Substitute Equations (14) and (16) into Equation (15a), we have. 


At t rci rri tfKi t 

/ (1 „ L )( !£i {T } _ + {T + (1 _ L) [K] {T > - {P .» dt = 0 

Q At At J At 1 At J At J 


After integration, it yields, 


,[C] At [C] At it 

("** * “ DO) tTj> = (-“ - -- DO) lij} + -- [?i] 


2 3 


(17) 


The constants in parentheses are weighting factors, which my be modified to 
improve the stability of the numerical solution process. Based on Donea’s 
suggestion [6], Equation (17) is modified to the following form. 



At [C] 

+ r DO) {T i> = c~ 


At At 

-- DU) {T-j } + ~ [Pi] 


( 18 ) 


Equation (18) is the equation of solution in time domain. 



The initial temperature of the blade is assumed to be uniform and same 
as the ambient temperature, 90 "F, namely, 

= {9 ° D} 

then, an iterative process is commenced to determine the temperature at the 
next time step t +At, namely (Tj}, by using Equation (IS). This process is 
repeated until the maximum difference between {Tj} and (Tjl is smaller than 
a predetermined small quantity. In this analysis, 500 second is chosen 
as A t and good results are obtained as shown in Figs. 13 and 14. 

Observing Figs. 13 and 14 and the computer results, numerical instability 
is experienced in the first hour, but thereafter, stability is restored and 
steady-state condition is reached. They indicate that the length of sunshine 
time required for the blade to reach the steady-state temperature is in the 
neighborhood of 4 hours. Obviously, this is possible. Thus, steady-state 
temperature distribution of the blade is used as the input of Phase II, finite 
element thermal -stress analysis. 


PHASE II. THERMAL-STRESS ANALYSIS 


A. PI ana-Strain Thermoelastic Problem 

From the previous discussions, it has been stated that the blade is a 
slender body, its length is large compared with its maximum cross-sectional 
dimension* The thermal input along the blade axis is close to uniform, 
namely, the blade temperature is independent of the axial coordinate. Under 
these conditions, the use of the concept of plain strain not only provides a 
good solution E7], but also greatly reduces the computational effort, because 
it decreases a three-dimensional problem to a two-dimensional one. 

If the orthogonal axes system e, n and € are used, the fundamental 
assumptions for a plane-strain thermoelastic problem are: 

- 0, - u ? U,ri), Mt, - (19} 

The corresponding strain components have the form: 

e S - Hi ~ e !n = 0 

e 55 “ e 55;(5, T i)» e nn ~ e nn^» ^ E ?n = ^0) 

The stress field is related to the strain field by the elastic constants of 
modulus of elasticity and Poisson's ratios. As stated before, wood is an 
orthotropic material, the modulus of elasticity of wood in longitudinal 
direction (parallel to grain) have in general, a significant difference with 
those in radial or tangential direction, some twenty times as large [9]. 

To account for the possible difference in radial and tangential directions, 
a 10% difference is again maintained in this computation. For plywood, an 
average value is taken. They are listed in Table 3. 


Table 3. 


Douglas Fir 

Birch 

Birch Plywood 

el 

1.95x10® psi 

2.10x10® psi 

1.38x10® psi 

e t 

0.098X10® psi 

\ 

0.105x10® 

0.77x10® 

Er 

0.133x10® psi 

0.164x10® 

0.164x10® 

g tr 

0.014x10® psi 

0.036x10 s 

0.076x10® 


The Poisson s ratios V and the thermal expansion coefficients, a, of 
wood are also different significantly in the major axis directions. They [9] 
are listed in Table 4. 


Douglas Fir 


Table 4. 
Birch 


Birch Plywood 


°L 

2.1xl0“®/°F 

2.1xl0®/°F 

14.57x10“® 

“r 

34.9X10"®/ D F 

39.52xl0-®/°F 

27.05x10"® 

«R 

25. 9x10"®/ °F 

30.38xl0"®/°F 

30.38x10"® 

V LR 

0.292 

0.426 

0.433 

V LT 

0.449 

0.451 

0.308 

V RT 

0.390 

0.697 

0.476 

V TR 

0.287 

0.447 

0.440 

V RL 

0.020 

0.033 

0.255 

V 

0.022 

0.023 

0.166 
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B. Formulation of Potential Energy 

In Phase I, the heat conduction portion of the problem, a governing 
partial differential equation needs to be established first, then, the finite 
element method in Variational approach is employed as a numerical procedure to 
obtain its solution. Now, for the problem of finding the thermal-stresses due 
to a temperature distribution, no controlling differential equation can be 
found. Thus, a conceptually different technique must be used to determine 
directly the displacement field, {fl, without going through the intermediate 
step of establishing the differential equation. Once {f> is found, the stress 
field {a} can be determined accordingly. 

First, with an assumed displacement field and under the influence of 
temperature rises at various locations, the total potential energy, P.E., may 
be expressed as 

P.E. = / (; [sJTtEl [s] - [eKE] C ctT] }dv (21) 

V 2 

where V represents the total volume of the blade. In a two-dimensional case, 
it may be replaced by the multiplication of thickness and cross-sectional area, 
i.e., t * A. 

Let v and w be the displacements of a point in ? and n directions, we may 
express [e] in terms of the displacement field {f} as follows: 

[a] {f} 


e 5 

E n 

, £8] = 

a 

a? 

0 

, if) = 

V 



0 

3 

an 


w 



a 

a 





an 





where 


(e> = 


Further, the strain components of a point, {e}, may be expressed in terms of 
nodal displacements {D}. For demonstration purpose, a triangular element as 
shown in Fig. 15 and the aforementioned shape functions, [N], are again used, 
the strain field, {e}, may be expressed as: 


or 


where 


{e} = [3] DU {01, 


{El - m {Dl 

CBl - [3] DU 


CH3 = 


Nj 0 Nj 0 n k 0 

0 Hi 0 N j 0 N k 


m - 



y j 


w* 


Vk 



( 22 ) 


!ii = “ Cai + bj e + Cj n] 

Nj 15 2k 5 * C J ^ 

1 

h = H tajc + b k ? + c k n] 

Note that A, a-f, aj and have been described in (lib) and (11c). 

Expanding {Dl to cover all the nodal displacements in the entire blade section 
and substituting Equation (22) into Equation (21) r , one obtains: 


P.E. = i [D]T £/[B]T[E][B]dvJ CD] - [D]T J[B]T[E3CaT3<lV (23) 


Equation (23) expresses the total potential energy of the blade section in 
finite element formulation. 


C» Finite Element Method— Minimum Potential Energy Approach 

Recall that the total potential energy, P.E f , is formed based on an 
assumed displacement field If), for each element. These fields are then 
expressed in terms of unknown nodal displacements by using shape functions. 
Although admissible fields are chosen, the compatibility conditions are 
therefore satisfied within the elements and at the nodes, yet the equilibrium 
conditions remain to be determined. 

It is known that among all admissible configurations of a conservative 
system, those that satisfy the equations of equilibrium make the potential 
energy stationary with respect to small variations of displacement. If the 
stationary condition is a minimum, the equilibrium state is stable. Thus, 
the equilibrium conditions may be sought through the minimization of the 
total potential energy P.E., namely. 


3D 


or, (Sj [B]T[E][B]dvl (DJ = | / [Bl T [E][aT3dv (25) 

1 V - V 

where v is the volume of an element. 

Equation (25) may be written in the following simplified form: 

DG [D] = CP 3 (26) 

where [K] is the total stiffness of the blade and [k] is the stiffness 
of each element, tp e ] is the thermal load produced at nodal points in a 

given temperature field and [E], Matrix of modulus of elasticity^ 
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RI 

DC] - E M 
1 

Oil = /CB]TCElCBldV 
£P] - S [p 6 ] 


(27) 


Cp e 3 - / [B3TCE]Cat]dv 

V 

Equation (26) is the finite element formulation of the blade section 
if isotropic material is used. But the wood composite is an ortho tropic 
material, the matrices of modulus of elasticity CE] in Equation (25), need 

t 

to be modified accordingly as to be described in the following sections. 


D. The Treatment of Orthotropic Property 

For an orthotropic material such as wood [8] and in the case of plain 
strain^the normal strain-stress relationships apply: 


also 


Thus 







On + 


1 





cf 0 

5 + -- v^n n 

E n 


(28) 


( 29 ) 


Note that v ?T] denotes the strain in the ? direction due to strain in 


the n direction. 
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Substitution of Equation (29) Into (28) produces: 

>L v 

S . i- (1 - v ?5 ' J K ) o e - -- (v 5 n + 5n ) ff n 

e n = ~ Kc + v nE V 55' °5 + i; 11 ‘ ° n 

Solving the above equations, one finds: 

^ « (aE^ + bE c %) / (ad - be) 

a n = (cE n % + dE n S n ) / (ad - be) 

_ r Y 
T ?n ~ b cn Sn 

, v 

where, a - 1 " v nC Sn 

b = v so + v a % 

c = * v n£ 

d - 1 - 
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(30) 


(31) 


or, the stress-strain relationships my be described In a matrix form as 


{cr> = [E m ] (e> 


(32) 


where 


[E m 3 


1 

ad - be 


% 
c E n 
0 


b E c 0 

d E n 0 

0 (ad - bc)G ?T1 


(33) 


From the Maxwell 


or. 


-Betti's Reciprocal 
cE ri = bE £ 
En b 


theorem, 


so the matrix of modulus of elasticity remains symmetrical for wood 
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Equation (33) gives the matrix CE m !l for plane strain in an orthotropic 
material, where c, n and ? are principal directions of orthotropy. 

The elasticity constants of wood, v, e, « and 8, have been described in 
Tables 3 and 4, In view of the blade construction, wood property axis T 
corresponds to axis ? in this section, R to n and L to €* 

E, Element Axes to Structural Axes 

Because of the blade curvature, the principal axes of material at one 
point in the blade section have in general, a different orientation with 
respect to the structural axes than those at another point. A typical 
orientation of laminated wood has been shown in Fig, 5 0 where s, n and K are 
the principal axes of material and local axes of the element, while y, z and x 
are the structural axes and the axes of the blade section. The transformation 
from material axes to structural axes is necessary for all elements in order 
to unify them in the framework of structural axes. This transformation process 
may be described as follows: 

T ms - Transformation from structural axes to material axes 


then, 

t E m^ = £Tu, s ] leg) (34) 

where subscript m and s represent material axes and structures axes, and 

cos 2 B sin 2 3 

sin 2 3 cos 2 0 


£ T ms3 - 


-2 sin 3 cos 


2 sin 0 cos 0 


sin 3 cos 3 
-sin 3 cos 3 
cos 2 3 -sin 2 3 


In addition, the stress transformation from material axes to structural 
axes is 

W = ^rns-l ^■ t W' 

A proper stresses train transformation may then be obtained by first stating 
the stress-strain relationships in material axes, namely, 

torf - ttml t e m } (361 

where [E m ] has been described in Equation (33). Substitution of Equation 
(34-) and (36) into (35) gives: 

{o s l = CT ms ] [E m l [T ms ] {e s } 
or 

{ts£ = [E s ] U 5 ) (37) 

where 

CE S ] = CT ras ]T [£„] CT ms a 

CE S 3 is the matrix of modulus of elasticity for orthotropic material, of which 
the material axes £ and n are oriented at p angle with respect to principal 
axes y and z, and material axis 5 coinsides with structural axis x. This 
matrix is used in place of [E] in Equation (25) for the determination of the 

stiffness of each element, CO, and thermal load matrix of each element, CP e 3. 

Then, proper combinations of DO and CP e 3 respectively for all elements render 
EK3 and [Pi as described in Equations (25) to (27). Finally, the equation 
(26) is solved and the displacements of all nodal points [0] are determined. 
When the displacement field is found, the stresses at each element can be 
computed accordingly. Thus, the determination of thermal-stress distribution 
on the wood blade is accomplished. 



F. Thermal Stresses 

Corresponding to the two steady-state Heat Conduction analyses made in 
Phase I, two thermal-stress analyses are run on the University computer. The 
stresses are all computed in material axes 5 and n as shown in Fig. 5. Only 
the stresses of elements in the surface layer are plotted. For example, 

Fig, 16 shows the stresses in g direction when the blade section is in 
horizontal position, and Fig. 17 shows the same type of stresses when the 
blade section is in vertical position. 

Observing Fig* 16 and 17, one finds that the high temperature sides of 
the blades, stresses of same magnitude are induced. This offers an indepen- 
dent check on the correctness of the computations. The sense of the stresses 
on the high temperature side is tensile, this is largely due to the special 
characteristics of the blade construction. As described before, in D-spar a 
thin layer of 1/8" Birch plywood wraps around 3" of lamintated Douglas Fir. 

The thermal expansion coefficient In tangential direction of Douglas Fir is 
about Z%% higher than that of Birch plywood. When the temperature rises, 
the thin layer of Birch plywood is 'stretched' by the expanding Douglas Fir. 
Thus, tensile stresses are produced on the exterior surface of the blade. 

This computational phenomenon agrees to the rationality of engineering. 
Furthermore, in view of the smooth transition from high stress area to low 
stress area; in Fig. 16 the stresses on the tail portion varies when the 
interior materials changed from Douglas Fir to honeycomb paper core; and the 
approximately symmetrical stress distribution on Fig. 17 for an approximately 
symmetrical blade geometry and thermal loadings; one can conclude, short of 
experimental verification, that the analysis offers a set of results with high 
believability. 



Conclusion Remarks 


A. Pre-run Program Testings 

For the purpose of increasing the reliability of the finite element 
analyses, both the heat conduction analysis computer programs and the thermal"* 
stress analysis computer programs have undergone several testings before using 
on the analysis of the blade. These testings consist of solving smaller but 
similar engineering problems with either theoretical solutions or solutions 
given by other authors, Nearly identical results have been obtained in all 
testings. 

B. Finer-Mesh Testing 

After the thermal-stress analysis has been performed on the blade, for a 
finite element model of 388 elements, the analysis is repeated on a finer-mesh 
model of 888 elements* For the finer-mesh model, the stresses in C direction 
of the surface layer elements are shown in Fig, 18. A comparison of Fig. 16 
to Fig. 18 shows that the stresses of the former are converging to the 
corresponding stresses of the later. This is considered to be a good 
convergence test for finite element analysis. 

C. Engineering Judgement 

The final results of steady-state temperature distribution and the 
thermal -stress distribution as shown in Figs. 6, 7, 16 and 17 are agreeable 
to engineering rationality as stated in section E of Phase I and section F 
of Phase II. 


D. The Cause of Cracking 

The allowable stresses in principal directions of Birch plywood [9] are 
listed in Table 5. A comparison reveals that the maximum blade thermal" 
stress, which occurs in C direction, is far lower than the allowable stress 
in the corresponding direction, Although the trancient thermal -stress 
analysis is not performed due to time limitation of the research project, yet 
it is believed that they should not exceed the allowable stress level for a 
prolonged period, since the sun is a mild heat source. 

Table 5 


Allowable Stresses for Birch Plywood 


Compressive 

* 

Tensile 

a L 

5,770 p s j 

13,640 p s -j 

°T 

3,370 

7,280 

°R 

970 

920 


E. Conclusion 

In conclusion, the cause of cracking of the blade is not due to thermal 
stresses induced by the sun. Investigation of other causes, such as poor 
manufacturing, needs to be made. In addition to the elimination of thermal 
stress as a possible reason of blade cracking, the research brings forth a 
method and a practical example of thermal-stress analysis for an engineering 
body of orthotropic materials, which will find many applications in the 
Engineering Community. 
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Figure 9, Temperature Distribution of Exterior Elements 
Blade Section in Horizontal Position 
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Figure 10. Temperature Distribution of Exterior Elements 
Blade Section in Vertical Position 
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Computer Programs and Output Examples 

(A) Grid Generation Program for Triangular Element A1 

(B) Bandwidth Reduction Program A6 

(C) Steady State Heat Conduction Program - Heat 1 A10 

(D) Transcient State Heat Conduction Program - Heat 2 A22 

(E) Changes on Grid Generation Program, Program (A), 

from Triangular Element to Rectangular Element A37 

(F) Changes on Bandwidth Reduction Program, Program (B), 

from Triangular Element to Rectangular Element A39 

(G) Transfer Program - Preparing Output Data from Programs 

(A) and (C) in the Form of Input Data of SAP IV A42 

(H) SAP IV Program for Thermal Stress Analysis when the Blade 

Section is in Horizontal Position - Rectangular Elements A45 

(I) Transformation of Stresses in Structural Axes of SAP IV 

to Stresses in Material Axes A65 

{ J ) Output Example 1: Temperature Distribution when the 

Blade Section Is in Horizontal Position A66 

(K) Output Example 2: Thermal Stresses of Elements Along 

Boundary When the Blade Section is in Horizontal Position A69 



original mm m 

OF POOR QUALITY 

c (<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 


C[ <> 
cj THIS PROGRAM GENERATES THE GRIDS FOR ANY GIVEN <> 
C j REGION <> 
cj TRAINGULAR FINITE ELEMENT <> 

cj <> 


cj <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 

c 

c 

DIMENSION N00(250Q,6),IRR(500) 

DIMENSION XPC500),YP(500) 3 XKG(9),YRG(9) 3 N(8) 3 NDN(8) 3 XC(21 3 21) 3 
1JT(74,4) } IC0MP(4,4) 3 NN(21,21) ,NNRB(74,4,21) ,XE(3000) ,YE(3G00) , 
1NE(3000) ,NRC4) ,LB(3) 3 TITLE(10) , IIRO(90) , IICQ(90) , 

2NOU(3000) 3 NOUK(3000) 3 YC(21 ,21) 

REAL N 

DATA NBW/O/ } NB/0/ ,NEL/0/ 

DATA ICOMP/-l,l, 1,-1, 1,-1, -1,1, 1,-1, -1,1, -1,1, 1,-1/ 

DATA IN/5/ , 10/6/ 

DATA IP l Ml 

C 

C FINITE ELEMENT GRID GENERATION PROGRAM 

C 

G<><><> <><><><><><><><><><><><><><> <><><><><><><><><><><><><><><><><><> 

C 

C TITLE=AN IDENTIFICATION STATEMENT 

C INRG=THE NUMBER OF REGIONS ORIGINALLY DEFINED 

C (CAN SPECIFY A MAX OF 20 REGIONS) 

C INBP=THE NUMBER OF BOUNDARY POINTS OR GLOBAL POINTS ORIGINALLY DEF 
C (MUST HAVE 8 IN EACH REGION) 

C IPCH=0 IF ELEMENT DATA IS NOT TO BE PUNCHED 

C IPCH=1 IF ELEMENT DATA IS TO BE PUNCHED 

C XP(I)=THE X COORDINATE OF EACH OF THE ORIGINAL BOUNDARY POINTS 

C YP(I)=THE Y COORDINATE OF EACH OF THE ORIGINAL BOUNDARY POINTS 

C NRG=THE REGION NUMBER 

C JT(NRG,J)=TKE REGION CONNECTIVITY DATA 

C NROWS=THE NUMBER OF ROWS OF NODES IN A REGION 
C NCOL=THE NUMBER OF COLUMNS OF NODES IN A REGION 

C (CAN SPECIFY A MAX OF 21 ROWS AND 21 COLUMNS IN A GIVEN REGION) 

C NDN(I)=THE GLOBAL NODE NUMBERS ORIGINALLY DEFINED IN A GIVEN REGIO 

C 
c 
c 

100 FORMAT ( 10 A4) 

101 FORMAT (3 13) 

102 FORMAT (8F6. 2) 

103 FORMAT C5I3,2A4) 

1113 F0RMAT(1X,5I10) 

104 FORMAT(////1X,10A4//1X,18HGLOBAL COORDINATES //IX, 

130HNUMBER X COORD Y COORD ) 

105 F0RMAT(2X, 13 3 7X 3 F7.4 3 5XjF7 . 4) 

106 FORMAT (// IX 3 17HC0NNECTIVITY DATA/ IX , 41HREGI0N SIDE 1 2 

1 3 4 ) 

107 FORMAT (2X, 13, 14X 3 4(I2 3 5X)) 



ORIGINAL Pft©E ES 
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108 F0RMAT(11I3) 

109 F0RMAT(///1X,12H*** REGION ,13,6k ***//10X,I2,5H R0WS,10X,I2 

1j8H COLUMNS // 1 OX , 2 1HB0UND ARY NODE NUMBERS, 10X, 815) 

110 F0RMAT(//1X, 19HREGI0N NODE NUMBERS/) 

111 F0RMAT(1X,2QI5) 

112 F0RMAT(//5X,17HNEL NODE NUMBERS , 9X , 4HX C 1 ) , 6X , 4HY ( 1 ) , 8X , 4HX ( 2 ) , 8X , 

14HY(2),8X,4HX(3),8X,4HY(3) ) 

113 FORMAT ( IX , 5 15 , 3X , 6F 12 . 4) 

114 FORMAT (41 4 , 6F14 . 5 ) 

115 FORMATC/// IX, 21HBANDWIDTH QUANTITY IS,I4,22H CALCULATED IN ELEMENT 
1,14) 

READ (IN, 100) TITLE 

C READ ( IN , 129 ) NKA , SEL , RMO , RA , RMI , SHR 

129 FORMAT(I3,5FlO .5) 

DO 400 IIM=1,1 

NEL=0 

NB=0 

READ(IN, 101) I NRG, INBP, I PCH 
WRITE (10, 116) INRG , INBP , IPCH 

116 F0RMAT(1X,3I3) 

RE AD (IN, 102) (XP (I), 1=1, INBP) 

WRITE(IO, 102) (XP(I), 1=1, INBP) 

READ (IN, 102) (YP (I), 1=1, INBP) 

WRITE (10, 102) (YP(I), 1=1, INBP) 

DO 1 1=1, INRG 

READCIN, 103) NRG, (JT(NRG,J),J=1,4) ,IRR(I) 

WRITE(I0, 117) NRG, (JT(NRG, J) ,J=1,4) , IRR(I) 

117 F0RMAT(//,1X,5I3,2A4) 

1 CONTINUE 

C WRITE (10, 104) TITLE 

IF(IIM.EQ.2)G0 TO 401 
WRITE (IP, 100)TITLE 

C WRITE ( IP , 129 ) NRA , SEL , RMO , RA , RMI , SHR 

401 WRITE (10, 105) (I,XP(I) ,YP(I) ,1=1, INBP) 

G%% WRITE (10, 106) 

DO 2 1=1, INRG 

G%% WRITE (10, 107) I, (JT(I, J) , J=l,4) 

2 CONTINUE 

C 


C* 

c 




C LOOP ON THE REGIONS IN ORDER TO GENERATE THE ELEMENTS 


.¥% <V #V ■*> <V 4k #% 4k 4 


JL. ki. kl- k.L.k£kk.l>klk kL JL k.1 


-I* JL Jt— 




DO 22 KK=1 , INRG 

READ(IN, 108) NRG, IIRO(KK) , IICO(KK) , (NDN(I) ,1=1,8) 
WRITE (10, 109) NRG, IIRO(KK) ,IICO(KK), (NDN(I) ,1=1,8) 
NR0WS=IIR0(KK) 

NC0If=IIC0(KK) 

C 

C GENERATION OF THE ELEMENT NODAL COORDINATES 
C 


A£ 


DO 3 1=1,8 


c 

c 

c 


c 


c 


II=NDN(I) 

XRG(I)=XP(II) 

3 YRGCI)-YP(II) 
XRG(9)=XRGC1) 
YRG(9)=YRG(1) 
TR-NROWS-1 
DETA-2 . O/TR 
TR=NCQL-1 
DSI=2. Q/TR 
DO 4 I=l,NROWS 
TR=I - 1 

ETA«1.0-TR*DETA 
DO 4 J=l,NGOL 
TR=J-1 


pm r 

OF PGQS QUALfr/ 


SI=-1.0+TR*DSI 

N(l)=- . 25* C 1 . 0~SI)*(1 • 0-ETA)* (SI+ETA+1 . 0) 
N ( 2)« . 5 0* (X . 0 -S I**2 )* C 1 . 0 -ETA) 
N(3)».2S*C1.0+SI)*C1.0-ETA)*(SI -ETA-1.0) 
N C4)= . 5 0* ( 1 ■ O+SI)* ( 1 . 0 -ETA**2 ) 
NCS)«.2S*(1.0+SI)*C1.0+ETA)*CSI+ETA-1.0) 
N(6)=.50*(l . 0 -SI**2)*(1 . O+ETA) 
N(7)=.25*C1.0-SI)*C1.0+ETA)*(ETA-SI-1.0) 
N (8)= . 5 0* C 1 • O-SI )* C 1 . 0-ETA**2 ) 

XC(I,J)=0.0 
YC(X, J)=O.Q 
DO 4 K=l,8 

XC (I , J)=XC (I , J)+XRG(K)*N CK) 
YC(I,J)«YC(I,J)+YRG(K)*N(K) 

4 CONTINUE 


GENERATION OF THE REGION NODE NUMBERS 


KN1=1 
KS1-1 
KN2=NRGWS 
KS2=NGOL 
DO 11 1=1,4 
NRT=JT(NRG, I) 

IFCNRT.EQ. O.OR.NRT.GT.NRG) GO TO 11 
DO 5 J=l,4 

5 IF(JT(NRT,J) .EQ.NRG) NRTS=J 
K=NCOL 

IFCI.EQ.2.0R.I.EQ.4) K=NROWS 
JL=1 

JK=ICOMP C I > NRTS ) 

IF(JK.EQ.-l) JL=K 
DO 10 J=1,K 
GO TO (6 , 7 , 8 , 9) , I 

6 NN(NROWS,J)=NNRB(NRT,NRTS,JL) 

WRITE (6, 7000 )NN(NROWS,J) 
KN2=NR0WS-1 
GO TO 10 

7 NN C J , NCOL)=NNRB (NRT , NRTS , JL) 

WRITE (6 , 700 1 ) NN( J , NCOL) ' 


AS 


n a a o n a cj a a 


owmmi mm m 

OF POOR QUAUW 


7000 

FORMAT f// ’ 
FORMAT (// ’ 
FORMATf//* 

STATION 

2 * 

, 19/) 

7001 

STATION 

1 1 

, 19/) 

7002 

STATION 

3 1 

, 19/) 

7003 

FORMAT f//' 

STATION 

4 ' 

, 19/) 


KS2-NC0L-1 
GO TO 10 

8 NN f 1 , J) =NNRB CNRT , NRTS , JL) 

C WRITE C6 , 7 004 ) NRT , NRTS , JL 

7004 FORMATf//’ NRT=',I9/' NRTS= , 3 I9/’ JL=’,19) 

C WRITE (6, 7002)NNfl, J) 

KN1=2 
GO TO 10 

9 NNfJ, 1)=NNRBCNRT,NRTS, JL) 

C WRITE(6, 70Q3)NNfJ, 1) 

KS 1=2 

10 JL=JL+JK 

11 CONTINUE 

IF(KN1.GT.KN2) GO TO 16 
IFfKSl.GT.KS2) GO TO 16 
DO 12 I=KN1,KN2 
DO 12 J=KS1,KS2 
NB=NB+1 

C WRITE f 10, 777)NB 

12 NNfl, J)=NB 

777 FORMATf////' NB = *,I6) 

STORAGE OF THE BOUNDARY NODE NUMBERS 

DO 13 I=1,NCQL 

NNRB (NRG , 1, I)=NN (NROWS , I ) 

13 NNRB fNRG ,3,1 )=NN (1,1) 

DO 14 1=1, NROWS 

NNRB fNRG ,2,1) =NN f I , NCOL) 

14 NNRBfNRG,4,I)=NNfI,l) 

OUTPUT OF THE REGION NODE NUMBERS 

WHITE f 10, 110) 

DO 15 1=1, NROWS 

WRITE flO, 111) fNNfl , J) , J=1 ,NCOL) 

15 GONTINUE 

DIVISION INTO TRIANGULAR ELEMENTS 

16 CONTINUE 
WRITE flO, 112) 

K=1 

DO 17 1=1, NROWS 
DO 17 J=1,NC0L 
KE(K)=XCfI, J) 

YEfK)=YGfI, J) 

NEfK)=NNfI,J) 

17 K=K+1 
L=NR0WS~1 


A4 


ORMAL PAST 0 
OF POOR QUALITY 

DO 21 1=1, L 
DO 21 J=2 , NCOL 

DIAG1=SQRT C (XC Cl , J) -XC CX+1 , J- 1 ) )**2+(YC ( I , J) -YG ( 1+1 , J- 1) )**2) 

DIAG2=SQRT( CXC C 1+1 , J ) ~XC ( I , J- 1 ) ) **2+ (YC(I+1, J) ~YC(I , J-l) )**2) 

NRC1)-NCOL*ItJ-1 

NR C 2 ) =NCOL* 1+ J 

NRC3)«NCOL*(I-l)+J 

NRC4)=NCOL*(I-l)+J-l 

DO 21 IJ=1 ,2 

NEL=NEL+1 

IF C (D1AG1/DI AG2) , GT . 1 . 02 ) GO TO 18 
J1=NR( 1) 

J2=NR(IJ+1) 

J3=NR(IJ+2) 

GO TO 19 

18 J1=NRCIJ) 

J2=NR(IJ+1) 

J3=NR(4) 

19 LBC13=IABS(NECJ1)-NECJ2))+1 
LB C2) =IABS CNE ( J2 ) ~NE C J3 ) )+l 
-LB C3)=IABS (NE ( J1 ) -NE C J3) )+l 
DO 20 IK=1 , 3 

IF(LB(IK) .LE.NBW) GO TO 20 
NBW=LBCIK) 

NELBW=NEL 

20 CONTINUE 

IF( KK . NE . 46 ) GO TO 888 
C WRITE C6, 666) J3 , NE(J3) 

666 FORMAT (///// ' J3 = ’,19,' NE(J3) - * ff 19 ) 

888 CONTINUE 

WRITE CIO, 113) KK,NEL,NE(J1) ,NE(J2) ,NE(J3) ,XE(J1) , YE(J1) ,XE(J2) , 

1 YE(J2),XECJ3),YECJ3) 

IFCIPCH.EQ.O) GO TO 21 

WRITE(IP,114) NEL,NE(J1),NE(J2) ,NECJ3) ,XE CJ1) ,YE( Jl) ,XE{J2) ,YEC*I2) 
1,XECJ3),YE{J3) 

NOOCNEL, 1)=KK 
NOOCNEL, 2 )=NEL 
NOOCNEL, 3)=NE(J1) 

NOOCNEL, 4)=NE(J2) 

NOOCNEL, 5 )=NECJ3) 

NOOCNEL, 6)=IRRCKK) 

21 CONTINUE 

22 CONTINUE 

WRITE CIO, 115) NBW,NELBW 
DO 95 IB=1 , 1000 
95 NOU(IB)=0 

C%% WRITE CIO, 35) 

3 5 FORMAT C / / / , I OX , ' **** B OUND ARY NODE NUMBERS **** * ) 

40 FORMAT C22I 3) 

DO 30 1=1 , INRG 

C%% WR ITE C 1 0 , 34 ) I , I IRO C I ) , 1 1 CO C I ) 

34 FORMAT C / 3X , T REGION ' , 1 2 , 3X , 1 2 , IX , ’ ROWS ' , 3X , 1 2 , IX , ' COLUMNS ' ) 

DO 30 11=1,4 
NRT=JTCI,II) 


AS 


ORffCJSNAL PA©g £! 
OF POOR QUALiW 


IF(NRT.GT.0)G0 TO 30 
K=IICO(I) 

IF(II.EQ,2.0R,II.£Q.4)K=IIR0(I) 

G %% WRITE (10 ,33) II, CNNRB(I,II,m),III=l,K) 

33 FORMAT (IX, 'SIDE' , 12 , 1 , 3 %., 1514) 

DO 90 I 11=1, K 
NII=NNRB(I , XI, III) 

90 N0U(NII)-1 
30 CONTINUE 

IKK=1 

DO 91 IB=1, 700 
IF(NOU(IB).EQ.O)GO TO 91 
NOUK(IKK)=IB 
IKK=IKK+1 

91 CONTINUE 
IKK=IKK“1 
WRITE (IP, 94) IKK 

94 FORMAT(/50X,I4,46X) 

WRITE (IP ,92) (NOUKCI ) , 1=1 , IKK) 

C%% WRITE (10,93) (NOUK( I ) , 1=1 , IKK) 

92 FORMAT (20 15) 

93 FORMAT (15 14) 

DO 333 1=1, NEL 

WRITE (IF, 575) ( NOO(I, J) , J=l,6) 

WRITE (15, 575 )( NOO(I, J) , J=l,6) 

C%% WRITE(I0,575)( N00(I,J) ,J=1,6) 

333 CONTINUE 

575 FORMAT (5 15, A4) 

576 FORMAT(1QX,5I10 , 10X,A10) 

400 CONTINUE 

STOP 

END 

//* 

//“ ==== == ===== == 

//* 

//* BANDWIDTH REDUCTION PROGRAM 
//* 

//* ■ ===== ■■■ =*=— = ■■ 

//* 

C 

COMMON/ AAA/ NOO ( 2 5 0 0 , 6 ) , NTT ( S 0 ) 

COMMON/BBB/ X (3 , 2000) , Y(3 , 2000 ) ,NB0 (1OQ0 ) , TITLE ( 10 ) 
CQMMON/CCC/ NODES , LMENTS , JT( 12000) ,HEMJT(24000) , JMEM(3000) , 
$ JNT(3000) , IDIFFjMINMAX 
DATA IN/ 11/, 10/ 6/, IP/ 12/ 

C 

READ (IN, 100)TITLE 
WRITE f IP, 100)TITLE 
100 FORMAT (10A4) 

122 FORMAT ( I 3 , 5F1Q , 5 ) 

J =1 
JJ=0 

150 READ (IN, 17) (JT(3000*(I-1)+J) , 1=1 ,3) , (X(I, J) ,Y(I , J) ,1=1,3) 
JT(9000-hJ)=0 


0 ! 

C! 

C! 

Cl 

: C! 

C! 
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OF FOOK 


NODES=MAXO(JT(J) , JT(3000+J) , JT(6000+J) , JJ) 

JJ-NODES 

IF(JT(3000+J).EQ,0)GO TO 152 
J=J+1 
GO TO 150 
152 LMENTS=J~1 
NODES=JJ 

12 FORMAT C/ / , 5X, 15HNUMBER OF NODES ,14, 15X, 

1 T NUMBER OF ELEMENTS * , 14 , f f ) 

WRITE (10 j 105) TITLE 
105 FORMAT (//, 1X,20A4,//) 

WRITE (10 , 12) NODES , LMENTS 
WRITE (10, 13) 

13 FORMAT ( // , 3X , 1SHNEL JT1 JT2 JT3, 10X,4HX(1) ,8X, 

$4HY(1) ,8X,4HX(2) ,8X,4HY(2) ,8X,4HX(3) ,8X,4HY(3)) 

17 FORMAT (4X , 314 , 6F14 . 5 ) 

DO 10 J=l, LMENTS 

C WRITE (10, II) J, (JT(3000*(I“1)+J) , 1=1,3) , (X(I , J) ,Y(I , J) ,1=1,3) 

11 F0RMAT(1X,4I5 ,3X,6F12.4) 

10 CONTINUE 

READ ( IN , 3 0 0 ) IBO 
300 FORMAT(50X, I4,46X) 

READ (IN, 301) (NBO (I) , 1=1 , IBO) 

WRITE (10, 302) 

DO 333 1=1, LMENTS 

READ (IN, 444) (NOO(I, J) , J=l,6) 

333 GONTINUE 

C WRITE (10, 303) (NBO (I), 1=1, IBO) 

302 FORMAT (//,' *** BOUNDARY NUMBERS ***') 

CALL SETUP 
NTBAN=IDIFF+1 
WRITE (10, 202)NTBAN 

202 FORMAT ( // , 2X , 25HTHE ORIGINAL BANDWIDTH IS, 14,//) 

CALL OPTNUM 

IF ( ID IFF . LE . MINMAX) GO TO 115 

MIBAN=MINMAX+1 

WRITE (10, 198) 

198 FORMAT ( 1H1/ / / , IX , 5 (3H0LD , 3X , 3HNEW , 7X) , / , IX , 

§5 (4HN0DE , 2X , 4HN0DE , 6X) ) 

WRITE (10, 200) (J,JNT(J),J=1, NODES) 

WRITE (7, 200) (J, JNT(J) ,J=1, NODES) 

200 FORMAT (5 (15 , IX, 15 ,5X) ) 

WRITE (10, 13) 

DO 180 J=l, LMENTS 
DO 555 11=1,3 
I I 1=11+2 

555 NOO(J, III)=JNT(JT(3000*(II-1)+J) ) 

C 

180 CONTINUE 

C WRITE(I0,11)J,(JNT(JT(3000*(I-1)+J)),I=1,3), 

.C $ (X(I , J) , Y(I, J) , 1=1,3) 

WRITE (10, 201) MIBAN 

201 FORMAT ( / / , 2X j 1 THE NEW BANDWIDTH IS', 13,/) 

WRITE ( I P , 3 10 ) NODES , LMENTS , MI BAN 


A 7 


ORIGINAL PASS m 
OF POOR QUALITY 


DO 181 J“1 , LMENTS 

181 WRITECIF, 18) J, (JNTCJTC3000*(I-1)+J)) , 1=1,3) , CX(I, J) , 

$Y(I,J) ,1=1,3) 

18 F0RMAT(4I5 , 6F10 .4) 

WRITE CIO, 302) 

WRITE (IP, 300) IBO 

WRITE (IP ,301) (JNT(NBO(I) ) , 1=1 , IBO) 

G WRITE(IO,303) (JNT(NBO(I)) ,1=1 , IBO) 

303 FORMAT (15 14) 

GO TO 117 
115 CONTINUE 

WRITE (10, 101) 

117 CONTINUE 

101 FORMAT (///,2X, 'THE ORIGINAL BANDWIDTH IS MININMUM* ) 

WRITE (IP , 3 1 0 ) NODES , LMENTS , NTBAN 
310 FORMAT (3 15 ) 

DO 190 J=1 , LMENTS 

190 WRITE(IP, 18) J, (JT(3000*(I-1)+J) ,1=1,3), 

1(X(I,J),Y(I,J), 1=1,3) 

WRITE (IP, 300) IBO 

WRITE (IP ,301) (NBO ( I ) , 1=1 , IBO ) 

301 F0RMAT(20I5) 

C— - C! 

444 F0RMAT(5I5,A4) 

LLL = N00(l,6) 

C 

LOL=0 

MMM=0 

NNN=1 

DO 666 I 1=1, LMENTS 

IF ( LLL„EQ,NOO(II j 6) )GO TO 1888 

NNN=NNN“ 1 

L0L=L0L+NNN 

MMM=MMMt1 

NTT(MMM)=LOL 

C WRITE(IO,888) NNN,NTT{MMM) 

LLL=NOO(II,6) 

NNN=1 

1888 CONTINUE 

C WRITE (10, 777 )NNN, (NOO(II, J) , J=l,6) , (X(I,II) ,Y(I,II) ,1=1,3) 
IF (II. NE . LMENTS) GO TO 1889 
LOL=LOL+NNN 
MMM=MMM+1 
NTT(MMM)=LOL 

C WRITE(IO, 888) NNN,NTT(MMM) 

1889 NNN=NNN+1 
666 CONTINUE 

C 

777 F0RMaT(4X,6I6,4X,A8,6F12.5) 

888 FORMAT (//’ NUMBER OF ELEMENTS IN THIS REGION = ! ,2110/) 
400 CONTINUE 
C 

C — - 


»jm£ WM 
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Vi 


STOP 

END 

SUBROUTINE OPTNUM 
C 

DIMENSION NEWJT(3000) , J0INTC3000) 

COMMON/ AAA/ N00 ( 15 00 , 6 ) , NTT ( 80 ) 

COMMON/ BBB/ X(3 ,2000) ,Y(3 ,2000) ,NBO(100n) ,TITLE(10) 
COMMON/CCC/ NODES, LMENTS, JT( 12000) , MEM JT ( 24 0 0 0 ) , JME M (300 0 ) , 
$ JNT(3000) , IDIFF.MINMAX 
MINMAX— IDIFF 
DO 60 IK=1, NODES 
DO 20 J =1, NODES 
JOINT (J)=0 
20 NEWJT(J)=0 

MAX=0 
1=1 

NEWJT(1)=IK 

JOINT(IK)=l 

K=1 

30 K4-JMEM (NEW JT ( I ) ) 

IF(K4.EQ.0)G0 TO 45 
JSUB- (NE W JT ( I ) - 1 ) *8 
DO 40 JJ=1 ,K4 
K5=MEM JT ( JSUB+JJ) 

IF(J0INT(K5).GT.0)G0 TO 40 
K=K+1 

NEWJT(K)-K5 
JOINT (K5)=K 
NDIFF=IABS(I~K) 

IF (NDIFF.GE. MINMAX) GO TO 60 
IF (NDIFF . GT . MAX) MAX=NDIFF 
40 CONTINUE 

IF (K.EQ. NODES) GO TO 50 
45 1=1+1 

GO TO 30 

50 MINMAX=MAX 

DO 55 J=l, NODES 
55 JNT(J)=JOINT(J) 

60 CONTINUE 

RETURN 
END 

SUBROUTINE SETUP 

COMMON/ AAA/ NOO(15O0, 6) ,NTT(BG) 

COMMON/BBB/ X(3 , 2000) , Y( 3 , 2000) ,NB0 ( 1000) , TITLE ( 10 ) 
COMMON/CCC/ NODES, LMENTS, JT( 12000) ,MEMJT(24G00) ,JMEM(3000) , 
$JNT(30Q0) , IDIFF, MINMAX 
IDIFF=0 

DO 15 J=l, NODES 
15 JMEM(J)=0 

DO 60 J=l, LMENTS 
DO 50 1=1,4 
JNTI=JT(3000*(I-1)+J) 

IF(JNTI.EQ.O) GO TO 60 
JSUB=(JNTI-1)*8 
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DO 40 11=1,4 
IF(II.EQ.I) GO TO 40 
JJT=JT(3000*CII-1)+J) 

IF(JJT.EQ.O) GO TO 50 
MEM 1=JMEM C JNTI ) 

IF(MEMl.EQ.O) GO TO 30 
DO 20 I I 1=1, MEM 1 

IF(MEMJT(JSUB+III).EQ.JJT) GO TO 40 
20 CONTINUE 

30 JMEMC JNTI )=JMEM ( JNTI )+l 

HEMJTC JSUB+JMEMC JNTI ) )=JJT 

IFCIABS ( JNTI- JJT).GT.IDIFF)IDIFF=IABS( JNTI -JJT) 
40 CONTINUE 
50 CONTINUE 
60 CONTINUE 
RETURN 
END 


C 

C 

C 

C 

C 

C&1B 
C 
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c <> 

REAL KXX,KXXX 

DIMENSION TITLE (10) ,NT1(11) ,NT2(80) 

DIMENSION NS (3) ,T(3) 

DIMENS ION IS IDE ( 2 ) , I CON ( 2 ) , PH I ( 3 ) 

DIMENSION CFV(6),N0DNUM(6) 

C 

INTEGER NP, IER,N 

REAL*8 2(482) ,R(482) ,AREAE(752) ,X(3) ,Y(3) ,DDT,TIME ,TDT 
2 ,X1 ,X2,X3 , Yl, Y2,Y3 ,COF,THSTF(482,35) ,THCAP(4S2, 35) 
COMMON/CONVEN/H ( 7 5 2 ) , I S IDEH (752,2) 

C 

COMMON / SOLID / Z , R , ARE AE , THSTF 
1 ,THCAP,NPT(752,3) 

REAL*8 ZZ (482 , 35 ) , WK(8694) ,F,EF 
COMMON / SOLVEQ / F(482, 1) ,EF(752) 

REAL *8 KZETA , KETA , RHQCP , BETA 

COMMON /CONT/ KZETA(752),KETA(752),RH0CP(752),BETA(752) 

C 

DATA IN/ 12/, 10/6/ 

C DATA NT1/343, 347, 351,355,359,363,367, 371, 375,379, 383/ 


C 

C DEFINITION OF THE CONTROL PARAMETERS 

C AREAE ( NE ) : ARRAY CONTAINS AREA OF EACH ELEMENTS 

C R ( NE*3): ARRAY CONTAINS THE X-CODRDINATES OF EACH ELEMENT 

C Z ( NE*3) : ARRAY CONTAINS THE Y-COORDINATES OF EACH ELEMENT 

C T ( 3 ) : ARRAY CONTAINS THE TEMP. OF THE ELEMENT'S THREE 


A to 
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C NODES. 

C TOLD ( NE*3): ARRAY COPIES THE T(3) ARRAY FOR EACH ELEMENT 

C NE = NUMBER OF ELEMENTS 

C NP - NUMBER OF GLOBAL TEMPERATURES 

C NBW = BAND WIDTH 

C H CONVECTION COEFFICIENT 

C TINF = AMBIENT AIR TEMP. 

C 

c 

c 

c 

C INPUT OF THE TITLE CARD AND THE CONTROL PARAMETERS 


TINF =10.00 
REWIND 12 


READ (IN, 3) 

TITLE 

WRITE (10,4 ) 

TITLE 

READ(IN,2) 

NP , NE , NBW 

NBW=NBW+17 


WRITE (10, 1102) 

NP,NE,NBW 


1 F0RMAT(4I5 ,6F10.4) 

2 F0RMATC3I5) 

3 FORMAT (10A4) 

4 F0RMAT(1H1,///1X,10A4) 
1102 FORMAT (IX, 316) 


OUTPUT OF TITLE AND DATA HEADINGS 


44 FORMAT ( 1H1 , / / / / IX , 20 A4/ / IX , 5HKXX =,F7 . 1, 10X,5HKYY = ,F7.1/1X, 
121HCONVECTION COEFF(l) =,F7. 1,3X,21HC0NVECTI0N C0EFF(2) =,F7.1/, 
2 IX, 15HFLUID TEMP(l) = , F7 . 1 , 3X , 15HFLUID TEMP(2) = ,F7 . 1// , IX, 
317HNEL NODE NUMBER , 6X , 4HX ( 1 ) ,6X,4HY(1) ,6X,4HX(2) ,6X,4HY(2) , 
46X,4HX(3) ,6X,4KY(3) ) 

A .U ^.JU 

ASSEMBLYING OF THE GLOBAL STIFFNESS MATRIX AND GLOBAL FORCE MATRIX 


INPUT AND ECHO PRINT OF ELEMENT DATA 


WRITE (6, 604) 

604 FORMAT (1H1, /////) 

WRITE (6 ,304) 

304 FORMAT (' T(l) T(2) T(3) 

1 ISIDE(l) ISIDEC2) '/) 

1061 F0RMATC6F10 . 3 , 2110) 

11 FORMAT ( * ' ,3F10.3,2I10) 

WRITE (6, 1107) 

WRITE (6 ,1108) 
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C NE : TOTAL NUMBER OF ELEMENTS 

C 

C. . IN: 12 MAIN 

T(1)=10, 

T(2)=10. 

T(3)=10. 

DO 180 1=1, NE 
H(I)=0.001D0 
ISIDEHCI , 1)=0 
ISIDEH(I ,2)=0 
180 CONTINUE 

DO 178 JJ=1,90 
N=2*JJ 

ISIDEH(N,1)=2 
ISIDEH(N,2)=0 
H(N)=5.D0 
178 CONTINUE 
C 

C BIRCH LOOP 

C 

DO 1001 J=l,280 
KZBTA(J)=. 0924 
KETA(J)~. 0836 

1001 CONTINUE 
C 

C DOUGLAS FIR LOOP 

C 

DO 1002 J=2S1, 302 
KZETA(J)=.0754 
KETA(J)=. 06820 

1002 CONTINUE 

DO 1010 J=35 1,638 
KZETA(J)=.0754 
KETA(J)=.06S2 
1010 CONTINUE 

DO 1015 J=687,752 
KZETA(J)=. 0754 
KETA(J)=. 0682 
1015 CONTINUE 
C 

C PAPER LOOP 

C 

DO 1020 J=303 ,350 
KZETA(J)=. 011161 
KETACJ)=. 011161 
1020 CONTINUE 

DO 1111 J =639,686 
KZETA(J)=. 011161 
KETA(J)=. 011161 
1111 CONTINUE 
C 

C LOOP ON THE REGIONS 

C 


LOOP 


a a 
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DO 7 KK=1,NE 

IFCKK.GEa.AND.KK.LE.22) BETA(KK)=-.2618 
IFCKK.GE.23.AND.KK.LE. 34) BETA CKK)=- . 15707 
IFCKK.GE.35.AND.KK.LE. 58) BETACKK)=0.0 
IFCKK.GE.59.AND.KK.LE.66) BETACKK)=. 38394 
IF CKK.GE.67.AND.KK.LE.74) BETACKK)=. 17452 
IFCKK.GE.75.AND.KK.LE. 82) BETA (KK)“1. 22164 
IFCKK.GE.83.AND.KK.LE. 90) BETA CKK)=, 69808 
IFCKK.GE.91.AND.KK.LE. 98) BETA CKK) =2. 443 3 
IFCKK.GE.99.AND.KK.LE. 106) BETACKK)=1. 91972 
IFCKK.GE.107.AND.KK.LE. 114) BETA CKK)=2. 9668 
IF(KK. GE.115.AND.KK.LE. 122) BETACKK)=2. 7574 
IF CKK. GE . 123 . AND . KK. DE . 146) BETA CKK)=3 . 14136 
IF CKK . GE . 147 . AND . KK . DE . 15 8 ) BETA C KK) =3 . 2 9 843 
IF CKK . GE . 159 . AND . KK . DE . 202 ) BETA CKK)=3 .40314 
IF CKK . GE . 203 . AND . KK . LE . 214) BETACKK)=3 . 29843 
IF CKK. GE.215.AND.KK.LE. 230) BETA CKK)=3 . 14136 
IF CKK . GE . 23 1 . AND , KK . LE , 246 ) BETA CKK) =0 . 0 
IFCKK.GE.247.AND.KK.LE. 258) BETACKK)=- . 15707 
IFCKK.GE.259.AND.KK.LE. 326) BETA(KK)=~.2618 
IFCKK.GE.327.AND.KK.LE. 350) BETA (KK)=- . 15707 
IFCKK.GE. 351.AND.KK.LE. 430) BETACKK)=0.0 
IFCKK. GE.431.AND.KK.LE. 438) BETA CKK) =.38394 
IFCKK.GE. 439. AND. KK.LE. 446) BETA CKK) =. 17452 
IF CKK. GE.447.AND.KK.LE. 454) BETA CKK) =.38394 
IF CKK . GE . 455 . AND . KK . DE , 462 ) BETA (KK)= .17452 
IF CKK. GE . 463 . AND .KK. LE . 47 0) BETA CKK) =1 . 2216 
IFCKK.GE.471.AND.KK.LE. 478) BETACKK)=. 69808 
IFCKK. GE.479.AND.KK.LE. 466) BETACKK)=1 . 2216 
IFCKK.GE. 487. AND. KK.LE, 494) BETA(KK)=. 69806 
IFCKK.GE. 495. AND. KK.LE. 502) BETA CKK) =2. 4433 
IFCKK. GE.503.AND.KK.LE. 510) BETACKK)=1. 91972 
IFCKK.GE. 511. AND. KK.LE. 518) BETA CKK)=2. 4433 
IF CKK. GE.519.AND.KK.LE. 526) BETACKK)=1 .91972 
IFCKK. GE.527.AND.KK.LE. 534) BETA CKK) =2. 9 66 8 
IF CKK . GE . 5 35 . AND . KK . LE . 542 ) BETA CKK) =2.75 74 
IFCKK. GE.543.AND.KK.LE. 550) BETA CKK)=2. 9668 
IFCKK. GE.551.AND.KK.LE. 558) BETA(KK)=2.7574 
IFCKK. GE.559.AND.KK.LE. 638) BETA CKK) =3. 14136 
IFCKK.GE. 639. AND. KK.LE. 662) BETA CKK)=3 . 29843 
IFCKK.GE. 663. AND. KK.LE. 708) BETACKK)=3.40314 
IFCKK.GE. 709. AND. KK.LE. 752) BETACKK)=0.0 

READ (IN, 1) NEL,NS ,X(1) , YU) ,XC2) ,YC2) ,XC3) , YC3) 

WRITEC6 3 121) NEL j NS 3 XC1),YC1) 3 XC2) > YC2),XC3) j Y( 3) ,BETACKK) 
C WRITE C7, 122) NEL.NS.XCl) ,YC1) ,X(2) ,YC2) ,XC3) ,YC3) .BETA CKK) 

121 FOEMATC4I5, 7F1G . 4) 

122 FORMAT C4I4 , 7F8 . 4 ) 

DIMENSION IN FEET. . . 

X1=XC1) 

X2=XC2) 

X3=XC3) 

Y1=YC1) 
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Y2=Y(2) 

Y3=Y(3) 

C 

X ( 1 )= X(l) /12.0D0 
X ( 2 )= X(2) /12.0D0 
X C 3 )= X(3) /12.0D0 

Y C 1 )= YCl) /12.0D0 

Y C 2 )- Y(2) /12.0D0 

Y ( 3 )= Y(3) /12.0D0 
C 

T1=T(1) 

T2=T(2) 

T3=T(3) 

1107 FORMAT (1H1,///20X 

* T(l) T(2) TC3) ' 

* ISIDEC1) ISIDEC2) ') 

1108 FORMAT C \2QX 

^ l 

3 - - - - ~ ~ ^ ^ “ 

^ t t *\ 

” 5 J 

C WRITE (10, 11)1, ISIDEHCKK, 1) , ISIDEHCKK, 2) 

C 

C 

C 

C TOLD(NS(l) )=T(1) 

C T0LD(NS(2) )=T(2) 

C T0LD(NS(3))=T(3) 

C H : CONVECTION COEFFICIENT 

C WRITEC6,5152)KK,NStl) 3 NSC2) J NSC3),T0LD(NS(l)),T0I,D(NS(2)) J 

C 1 T0LD(NS(3)), H(KK) 3 (ISIDECJJ) , JJ=1 ,2) , 

C 1 X1,Y1,X2,Y2,X3,Y3,N0T(KK,1),N0T{KK,6) 

5152 FORMAT ( 1 ' ,415, 3F10.3,F5. 0,212, 6F8. 4,15, A8) 

C 

T(1)=T1 
T(2)=T2 
T(3)=T3 

ECHO THE INPUT 

WRITE (10 } 23) NEL,NS jX(l) ,Y(1) ,X(2) ,Y(2) ,X(3) ,Y(3) 

* , HH,RHO,CP,ISIDE,T 

23 FORMAT(1HO,I5,1X,3I5,2X,6C1X,F8.4),2X,2F5.0,F6.3,2X,2I3, 

* 2X 3 3F6 . 0) 

C CALCULATION OF THE CONDUCTION MATRIX 
C 

C NPT (NEj 3) : ARRAY STORES THE NODE NUMBERS FOR EACH ELEMENT 
NPT(KK, 1) = NS(1) 

NPT(KKj 2) = NS (2) 

NPT(KKj 3) = NS (3) 

R(NSCl)) = X(l) 

R(NS(2)) = X(2) 

R(NS(3)) = X(3) 

ZCNSC1)) =Y(1) 

Z(NS(2)) =Y(2) 

Z (NS (3) ) =Y(3) 

AREAE (KK)=(X(2)*Y(3)-X(2)*Y(1)+X(3)*Y(1)“X(3)*Y(2) 
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1 + X(l)*Y(2)-X(l)*Y(3))/2.0 

— -0! 0! 

- Cl 

WRITEC6 ,266)KK,X(1) , Y(l) ,X(2) ,Y(2) ,X(3) , Y(3) ,AREAE(KK) 

266 FORMAT ( 1H0 , 15 , 6 ( IX , F8 . 4 ) ,2X3014.7) 

WRITE(6,319)KK,(T(I),I=1,3) 

319 FORMATC' ELEMENT NO. ' ,I8,3F12.4) 

7 CONTINUE 

END OF THE DO LOOP Cl 

WRITE(6 ,264) (I , TOLD (I) , 1=1 ,NP) 

264 FORMATC 1H0 , 15 ,E14 , 5 , 2X , 15 , E 14 . 5 , 2X , 15 , E14 . 5 , 2X , 15 ,E14 . 5 

* ,2X,I5,E14.5,2X,I5,E14,5) 

ADD EXTERNALLY APPL. CONC. NODAL FORCES VECTOR TO EF VECTOR 
ID1=0 
INK=0 

202 READ (IN, 2^3) NODNUM , CFV 

203 F0RMAT(6I3,2X,6F10.5) 

ID=0 

DO 204 L=l,6 

IF( NODNUM (L).LE.O) GOTO 205 
ID=ID+1 

INODNU=NODNUM(L) 

204 EF ( INODNU)=CFV (L) 4-EF ( INODNU) 

GOTO 206 

205 INK=1 
IF(ID.EQ.O) GOTO 216 

206 IF(ID1.EQ. 1)G0T0 222 
WRITE (10, 208) 

208 FORMAT (1H1, ' NODE NUMBER 1 ,10X, ' GENERALIZED NODAL FORCES ’) 

222 WRITE (10, 207) (NODNUM CL) ,CFV(L) ,L=1, ID) 

207 FORMATC 1H0, 6(13, E14. 5 ,2X)) 

IF(INK.EQ-1)G0T0 216 

ID 1=1 
GOTO 202 
216 CONTINUE 

C 

JJJKNT=0 

C 

CALL FORMTK (NE,NP,NBW) 

CALL FORMTF (NE,NP,NBW) 

CALL AS SMB L (NP,NBW,ZZ) 

C WRITE(6,9999) ((ZZ(I , J) , J=1,NBW) , 1=1 ,NP) 

9999 F0RMAT(2X,7G11.5) 

C WRITEC6, 9988) (F(I,1), 1=1, NP) 

9988 F0RMAT(2X, Gil ,5) 

CALL LEQT1B (ZZ , NP , 17 , 1 7 , NP , F , 1 , NP , 0 , WK , IER ) 

JJJKNT=JJJKNT+1 

IF ( JJJKNT . GE . 1 ) GOTO 979 

GOTO 980 

979 CONTINUE 

980 CONTINUE 
WRITE (10, 666) 
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666 FORMAT ( f NODE NO. * ,7X, ‘TEMPERATURE IN DEGREE F',/) 

WRITECe,787)(I,FCI,l),I“l,NP) 

787 FORMAT ( 2X, 15 , 10X,F15 . 8) 

WRITEC7,777)CF(I,1)»I«1,NP) 

777 FORMAT(9F8.4) 

C 

C 



STOP 

END 

SUBROUTINE FORMTK (NE ,NP,NBW) 

•jV****5V*********;V5V**5V*****************iVi‘f*' 

. FORM BANDED THERMAL CONDUCTANCE MATRIX 
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IMPLICIT REAL*4(A-H,0-Z) 

REAL* 8 KZETA , KETA , RHOCP , BETA 

REAL* 8 BI , B J , BK , Cl , CJ , CK , COF , CONI , CON2 , COM3 

RE AL*6 RB AR , CON , LG , RCTL , ELKT ( 3 , 3 ) , TKSTF (482,35) 

REAL*8 ADO , All , A22 , BOO , B 1 1 , B22 , COO , THCAP (482 , 35 ) 

REAL*8 Z (482) ,R(482) ,AREAE(752) ,X(3) ,Y(3) , DDT, TIME ,ATIM,P 
COMMON/CONVEN/H(752) ,ISIDEH(752 , 2) 


COMMON / SOLID / Z , R , AREAE , THSTF 
1 ,THCAP,NPT(752, 3) 

REAL*8 ZZ(482,35) ,F,B ,FF,EF 
COMMON / SOLVEQ / F(482, 1) ,EF(752) 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RHOCP(752) ,BETA(752) 

C 

C 

c 

C INITIALIZE GLOBAL THERMAL STIFFNESS MATRIX 

JN=5 

10=6 

DO 120 J=l,NBtf 
DO 110 1=1, NP 
THSTF(I,J) = 0 .DO 
110 GONT1NUE 
120 CONTINUE 
C 
C 

130 NBAND = 0 

DO 899 J=1,NE 

C WRITE(6,788)KZETA(J) ,KETA(J) ,BETA(J) 

788 FORMAT ( 3X , 3F20 . 9 ) 

899 CONTINUE 
C 

C WRITE (6 ,7776) 

7776 FORMAT ( * CONDUCTION MATRIX ELEMENTS ARE * ,/) 

DO 230 N=1,NE 

C<> <><><><><> — — — <><><><><><><><><><><> 


II = N.PT(N, 1) 

JJ = NPT(N, 2) 

KK = NPT(N,3) 

BI = Z(JJ) -Z(KK) 
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BJ = Z(KK)-Z(II) 

BK = Z(II)-Z(JJ) 

Cl = R(KK)-R(JJ) 

CJ « R(II) -R(KK) 

CK « R(JJ)-R(II) 

X(l) »R(II) 

X(2) =R(JJ) 

X(3) =R(KK) 

Y(l) =Z(II) 

Y(2) =Z(JJ) 

Y(3) =Z(KK) 

C WRITE (6,555)X(1),X(2),X(3),Y(1),Y(2),Y(3) 

555 FORMAT ( 2X , 6F 1 2 . 6 ) 

ALPHA=BETACN) 

C CON1=(KZETA(N)*(COS (BETA(N) ) ) **2+KETA ( N ) * (S IN ( BETA (N ) ) )**2) 

C 1/(4.*AREAE (N)) 

C CON2=(KZETA (N) -KETA (N) )*SIN(2 . *BETA(N) ) / (4 . *AREAE (N) ) 

C CON3=(KZETA(N)* (SIN (BETA (N) )) **2+KETA (N)*( COS ( BETA (N) ) )**2) 

C 1/ (4 . *AREAE CN) ) 

C0N1=(KZETA(N)*(C0S (ALPHA) )**2+KETA(N)* (SIN (ALPHA) )**2) 

1/ (4. "AREAE (N) ) 

CON2=(KZETA(N) -KETA(N) )*S IN ( 2 . *ALPHA ) / (4 .*ASEAE (N) ) 
CON3=(KZETA (N)* (SIN (ALPHA) )**2+KETA (N) -’'(COS (ALPHA) )**2) 
1/(4. •- ARE AE (N)) 

WRITE (6 , 7 9 99 ) CONI , C0N2 , C0N3 , AREAE (N) 

7999 F0RMAT(2X,3F10. 3, 2X,F10 . 7) 

ELKT(1 , 1) = (BI**2)*CON1+(CI**2)*CON3+(BI*CI)*CON2 
ELKT(1,2) = (BI*BJ)*C0N1+(CI*CJ)*C0N3+(BI*CJ)*C0N2 
ELKT(1,3) = (BI*BK)*C0N1+(CI*CK)*C0N3+(BI*CK)*C0N2 
ELKT(2, 1) = (BI*BJ)*C0N1+(CJ*CI)*C0N3+(BJ*CI)*C0N2 
ELKT(2 ,2) = (BJ— 2)*CONl+(CJ**2)*CON3+(BJ*CJ)*CON2 
ELKT(2,3) = (BJ*BK)*CONl+(CJ*CK)*CON3+(BJ*CK)*CON2 
ELKT(3, 1) = (BE*BI)*C0NX+(CK*CI)*C0N3+(BR*CI)*C0N2 
ELKT(3,2) = (BE*BJ)*C0N1+(CK*CJ)*C0N3+(BR*CJ)*C0N2 
ELKT(3 ,3) « (BK**2) *CONl+ (CK—2) *C0N3+ (BK*CK) *C0N2 
WRITE(6,*)N,H(N) , ( ISIDEH(N, Jl) ,31=1,2) 

DO 10 IQ=1,2 

C WRITE (6, 7805 )N, ( ISIDEH(N, Jl) , Jl=l,2) 

7805 FORMAT ( 1 SIDE EXPOSED ’,3110) 

IF ( ISIDEH(N,IQ) . LE . 0 ) GO TO 240 
JQ = ISIDEH (N,IQ) 

C MITE (10, 12) JQ,N 

12 FORMAT (1H0, 1 INSIDE 1 , 13 , 3X, T OF ELEMENT' , 15) 

KQ = JQ +1 

IF ( JQ.EQ . 3) KQ=1 

LG = DSQRT((X(KQ)-X(JQ))**2+(Y(KQ)-Y(JQ))**2) 

HL =H(N)*LG/6. 

MITE (6,*) N,H(N) ,LG,HL 
IF( JQ.EQ . 1) GO TO 20 
IF (JQ.EQ. 2) GO TO 30 
IF (JQ.EQ. 3) GO TO 40 
GO TO 10 

20 ELKT(1, 1)= ELKT(1,1)+HL * 2. 
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ELKT(2,2)= ELKT(2,2)+HL * 2. 
ELKT(1 ,2)= ELKT(1,2)+HL 
ELKT(2,1)=ELKTC2,1)+HL 
GO TO 10 
C 

30 ELKT(2,2)= ELXT(2,2)+HL*2. 

ELKT(3,3)= ELKT(3,3)+HL*2. 
ELKT(2,3)= ELKT(2 ,3)+HL 
ELKT(3 , 2)=ELKT(3 ,2)+HL 
GO TO 10 
C 

40 ELKT(1,1)= ELKT(1,1)+HL * 2. 

ELKT(3 , 3)= ELKT(3,3)+HL * 2. 
ELKT(1,3)= ELKT(1,3)+HL 
ELKTC3, 1)=ELKT(3, 1)+HL 
10 CONTINUE 


WRITE CIO t 2000)Nj II, JJ,KK,XC1) ,Y(1) ,X(2) , Y(2) ,X(3) , Y(3) , CONDTY(N) , 
1AREAE(N) ,RBAR,C0N,ELKT(1, 1) ,ELKT(1,2) ,ELKT(1 ,3) ,ELKT{2,2) , 
2ELKT(2,3) ,ELKT(3,3) 

FORMAT C 1H0 f 415 , 8G1 1 . 4 , / 2 IX , 8G11 . 4) 

0 ! 
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240 


7543 


210 

220 

230 


5050 


CONTINUE 

WRITEC I0,266)N,XC1) ,Y(1) ,XC2) ,Y(2) ,X(3) ,Y(3) ,AREAE(N) ,NBANB 
FORMAT C 1H0 , 15 , 6 ( IX ,FB . 4) , 2X ,D14 . 7 , 2X , 15 ) 

WRITE (6,277)N,(( ELKTCI , J) , J=1 , 3) , 1=1, 3) 

FORMAT ( ' ’ ,15/ 3(3P25.6/)) 

DO 220 LL-1,3 
I=NPT(N ,LL) 

DO 210 MM=1,3 
J-NPT (N , MM) - T-M8 

THSTF C I , J) =THSTF (I, J)-fELKT(LL,MM) 

WRITE (6, 7543) I , J, II , JH,NPT(N, J) ,NPT(N,I) 

FORMAT C f ',618, ' ****** 1 ) 

IF (J.GT.NBAND) NBAND=J 
CONTINUE 
CONTINUE 
CONTINUE 

WRITE C6 ,5050) C (THSTF ( I , J) , J=1 ,NBW) ,1=15,20) 

FORMAT (7G1 1.5) 


WRITE (6,200) NBAND 
200 FORMAT (1H0,' NBAND= ',15) 

END OF PROGRAM 


RETURN 

END 

SUBROUTINE FORMTF (NE , NP , NB W ) 
FORM THERMAL FORCE VECTOR 
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REAL*S ELF(3) ,LG,RCTL,C0F S THCAP(482,35) ,THSTF(482,35) 
REAL*8 Z(482) ,R(482) ,AREAE(752) ,X(3) ,Y(3) ,DDT 
COMMON/ CONVEN/ H (752) , ISIDEH(752 , ?) 

COMMON / SOLID / Z , R , AREAE , THSTF 
1 , THCAP , NPT (752,3) 

REAL*8 ZZ (482 , 35 ) , F , EF , B , ELF1 , ELF2 , ELF3 , FF , HLL 

COMMON / SOLVE Q / F(482 , X) ,EF(752) 

REAL *8 KZETA , KETA , RHOCP , BETA 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RHOCP(752) ,BETA(752) 
REAL *8 QR(752) ,HL(752) 

DATA PI/3. 14159265D0/ 

10=6 

TINF = O.DO 


FOR EACH ELEMENT .FORM THERMAL FORCE VECTOR 

WRITE (6 ,9999) 

9999 FORMAT ( ' FORCE VECTORS ARE ' ,/) 

DO 110 1=1, NE 
EFCD = O.DO 
QRCD=o.o 
110 CONTINUE 

DO 1070 J=1 ,45 
N=2*J 

QR(N)=326 . 0 
1070 CONTINUE 

DO 230 N=1,NE 
ELF(1)=0.D0 
ELF(2)=0.D0 
ELF(3)=0 -DO 
II=NPT(N, 1) 

JJ=NPT(N,2) 

KK=NPT(N, 3) 

Cl = R(KK) -R(JJ) 

CJ = R(II)“R(KK) 

CK = R(JJ)-R(II) 

X(1)=R(II) 

X(2)=R(JJ) 

X(3)=R(KK) 

Y(1)=2(II) 

YC2)=Z(JJ) 

Y (3)=Z (KK) 

DO 10 IQ=1,2 

IF(ISIDEH(N,IQ).LE.O)GO TO 240 
JQ=ISIDEH (N, IQ) 

KQ=JQ+1 

IFCJQ.EQ.3) KQ=1 

LG=DSQRT( (X(KQ) ~X(JQ) )**2+ (Y (KQ) -Y(JQ) )**2) 

C WRITE (6, 999 )LG 
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999 FORMAT(2X,F10.S) 

HL(N)=QR (N)+H(N)*TINF 
IF C JQ.EQ. l)GO TO 20 
IF ( JQ.EQ. 2)G0 TO 30 
IF (JQ .EQ.3JG0 TO 40 
GO TO 240 

20 ELF(l) -HL (N)*LG* . 5+ELF ( 1 ) 

ELF(2) "HL (N ) *LG* . 5+ELF ( 2 ) 

GO TO 10 

30 ELF (2) =HL (N )*LG* . 5+ELF (2 ) 

ELF(3) =HL(N)*LG* . 5+ELF (3) 

GO TO 10 

40 ELF(l) =HL(N)*LG*.5+ELF(1) 

ELF (3) =HL(N)*LG*.5+ELF(3) 

10 CONTINUE 

C 

240 EFCH)- EF(II)+ELF(1) 

EF(JJ)=EF(JJ) +ELF(2) 

EF(KK)=EF(KK) +ELF(3) 

707 FORMAT ( r 1 ,2I3,5F15 .5) 

C IF(EFCII).EQ.O) GO TO 50 

C IF(EF(JJ) .EQ.O) GO TO 50 

C IFCBF(KK).EQ.O) GO TO 50 

C GOTO 230 

C50 WRITE(6, 1999)N,II, JJ,KK,EF(II) ,EF(JJ) ,EF(KK) 

1999 FORMATClHO ,416 ,2X,3F10.4) 

230 CONTINUE 

C WRITE CIO 5 2000 )(EF Cl), 1=1, NP) 

2000 FORMAT C2X, FI 0.4) 

C WRITE CIO, 1000) 

1000 FORMAT ( 1H0 , * LEAVING FORMTF IN TDHEAT . . . IB ' ) 

RETURN 
END 

SUBROUTINE ASSMBL (NP,NEW t ZZ) 

C ASSEMBLE MATRICES FOR RECURRENCE FORMULAS 

REAL *8 Z(482) ,R(482) ,AREAE(752) 

REAL*8 DDT , DT, TKSTF C482 , 35 ) ,THCAP (482 , 35 ) 

COMMON/CONVEN /H ( 7 5 2 ) , IS IDEH (752,2) 

C 

COMMON / SOLID / Z , R , ARE AE , THSTF 
1 , THCAP , NPT (752,3) 

REAL*8 ZZ(482 ,35) ,F,EF ,B ,FF(15 , 1) 

COMMON / SOLVEQ / F(482, 1) ,EF(752) 

REAL *8 KZETA , KETA , RHOC P , BETA 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RH0CP(752) ,BETA(752) 

DO 150 1=1, NP 
C 

C COEFFICIENT MATRIX 

C 

DO 110 J=1,NBW 
ZZCI.J) = THSTF (I, J) 

110 CONTINUE 

C WRITE (6 ,2000) (ZZ(I , J) , J=12,26) 


$20 
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2000 FORM AT C ( ,6E9.5) 


original ptm m 
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F(I,1> = EFCI) 

150 CONTINUE 

WRITE (6, 1313) (F(I,1) ,I=1,NP) 
313 FORMAT ( * 1 ,G11.5) 

1001 FORMAT C1H ,13,11611.5) 

END OF PROGRAM 

RETURN 

END 
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sV*****;V********************5V****************i^*A*******iVi';****sV**iW:5Wf* 

<> 

C&1B HEAT . . . 2 • <> 

C 

C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 

c <> 

REAL KXXjKXKX 

DIMENSION TITLE (10) ,NT1 (11) ,NT2(80) 

DIMENSION NS(3) ,T(3) 

DIMENSION ISIDEC2) , ICONC2) ,PHI(3) 

DIMENSION CFV(6) ,NODNUM(6) ,FI J(463 , 100) , TIME (100) 

C 

INTEGER NP,IER,N 

REAL*8 Z(483) ,R(483) ,AREAE(752) ,X(3) ,Y(3) , DDT, TIME, TDT 
2 ,X1,X2,X3 ,Y1, Y2 ,Y3 ,COF,THSTF(483, 18) ,THCAP(483, 18) ,TOLD(483) 
COMMON/CONVEN/H ( 75 2 ) , I S IDEH (752,2) 

C 

COMMON / SOLID / Z ,R,AREAE .THSTF 
1 ,THCAP,NPT(752, 3) ,TOLD 

REAL*8 ZZ (483, 18), WK(IOOOO), F,EF,MAX(500) 

COMMON / SOLVEQ / F(483) ,EF(752) ,ZZ 
REAL *8 KZETA ,KETA , RHOCP , BETA 

COMMON /CONT/ KZETA (75 2) ,KETA(752) ,RH0C?(7S2) } BSTA(752) 

C 

DATA IN/ 12/, 10/ 6/ 

C DATA NT1/343, 347, 351, 355, 359, 363, 367, 371,375, 379, 383/ 
C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


DEFINITION OF THE CONTROL PARAMETERS 


AREAEC NE ) 

R ( NE*3) 
Z ( NE*3) 
T ( 3 ) 


ARRAY CONTAINS AREA OF EACH ELEMENTS 
ARRAY CONTAINS THE X-COORDINATES OF EACH ELEMENT 
ARRAY CONTAINS THE Y-COORDINATES OF EACH ELEMENT 
ARRAY CONTAINS THE TEMP. OF THE ELEMENT'S THREE 
NODES. 

TOLD ( NE*3): ARRAY COPIES THE T(3) ARRAY FOR EACH ELEMENT 
NE » NUMBER OF ELEMENTS 

NP = NUMBER OF GLOBAL TEMPERATURES 

NEW = BAND WIDTH 
II = CONVECTION COEFFICIENT 
TINF = AMBIENT AIR TEMP. 


INPUT OF THE TITLE CARD AND THE CONTROL PARAMETERS 


C DELTA T: TIME INTERVAL FOR EACH ITERATION 
TDT=500 . 0 
DDT=TDT/ 3600. 

REWIND 12 

C 

READ (IN, 3) TITLE 

C WRITE (10,4 ) TITLE 

READ (IN, 2) NP , NE , NBW 
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WRITE(IO, 1102) NP,NE,NBW 

C 

1 FORMAT (41 5 , 6F10 . 4) . 

2 FORMA1T3I5) 

3 FORMAT (10A4) 

4 FORMAT(1H1,///1X,10A4) 

1102 FORM AT( IX, 316) 

C 

C 

c 

C OUTPUT OF TITLE AND DATA HEADINGS 

C - 

44 F0RMAT(1H1 ,////lX,20A4//lX,5HKXX = ,F7 . 1, 10X.5HKYY - ,F7.1/1X, 
121HC0NVECTI0N COEFF(l) = ,F7 . 1, 3X,21HCONVECTION C0EFF(2) =,F7.1/, 
2 IX , 1 5HFLUID TEMP ( 1 ) =,F7 . 1 , 3X, 15HFLUID TEMP (2) -,F7.1//,1X, 

3 1 7HNEL NODE NUMBER , 6X , 4HX ( 1 ) , 6X , 4H Y ( 1 ) , 6X , 4KX (2) ,6X,4HY(2) , 

46X,4HX(3) ,6X,4HY(3) ) 

0 .’..T. .U’j, .U 

C ASSEMBLYING OF THE GLOBAL STIFFNESS MATRIX AND GLOBAL FORGE MATRIX 
C 

C INPUT AND ECHO PRINT OF ELEMENT DATA 
G 

C 

C WRITE (6, 604) 

604 FORMATC1H1, /////) 

C WRITE (6 ,304) 

304 FORMATS T(l) T(2) T(3) 

1 ISIDE(l) ISIDE(2) '/) 

1061 F0RMAT(6F10. 3,2110) 

11 FORMAT ( * ’ ,3F10. 3,2110) 

C WRITE(6,1107) 

C WRITE (6, 1108) 

C 

G NE : TOTAL NUMBER OF ELEMENTS 

C 

C. .IN: 12 MAIN LOOP 

T(l)=10. 

T(2)=10 . 

T(3)=10 . 

DO 180 1=1, NE 
H(I)=0.0D0 
ISIDEH(I,1)=0 
ISIDEH(I,2)=0 
180 CONTINUE 

DO 178 JJ~1,90 
N=2*JJ 

ISIDEHCN, 1)=2 
ISIDEH(N,2)=0 
H(N)=5.D0 
178 CONTINUE 
C 

C BIRCH LOOP 

C 


CF9G1K&. ■ 
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DO 1001 J-1,280 
KZETA(J)=.0924 
KETA(J)=.0836 
RH0CP(J)=10. 94246 

1001 CONTINUE 
C 

C DOUGLAS FIR LOOP 

C 

DO 1002 J=281, 302 
KZETA(J)=.0754 
KETA(J)=. 0682 
RH0CP(J)=8.5488 

1002 CONTINUE 

DO 1010 J=35 1,638 
KZETA(J)=.0754 
KETA(J)=. 0682 
RHOCP ( J)=8 . 5488 
1010 CONTINUE 

DO 1015 J=687 3 752 
K2ETACJ)=.0754 
KETA(J)~.0682 
RHOCP (J)=8. 5488 
1015 CONTINUE 
C 

C PAPER LOOP 

C 

DO 1020 J=303,350 
KZETA(J)=. 011161 
KETA(J)=. 011161 
RHOCP (J) =8, 1357 
1020 CONTINUE 

DO 1111 J =639 j 686 
KZETA(J)=. 011161 
KETA(J)=. 011161 
RHOCP CJ) =8. 1357 
1111 CONTINUE 
C 

C LOOP ON THE REGIONS 

C 

DO 7 KK=1,NE 

IF(KK - GE . 1 . AND . KK . LE . 22) BETA (KK)=- . 2618 
IFCKK.GE.23.AND.KK.LE.34) BETA(KK)=-. 15707 
IFCKK.GE. 35. AND.KK.LE. 58) BETA(KK)=0.0 
IFCKK.GE.59.AND.KK.LE.66) BETACKK)=. 38394 
lFCKK.GE.67.ANf .KK.LE.74) BETA CKK)=. 17452 
IFCKK.GE.75.AND.KK.LE.82) BETA (KK) =1.22164 
IFCKK.GE.83.AND.KK.LE.90) BETA (KK)=. 69808 
IFCKK.GE.91.AND.KK.LE. 98) BETA (KK)=2 .4433 
IFCKK.GE.99 .AND.KK.LE. 106) BETA(KK)=1 . 91972 
IFCKK.GE. 107.AND.KK.LE. 114) BETA (KK)=2. 9668 
IFCKK.GE. 115. AND.KK.LE. 122) BETA(KK)=2. 7574 
IFCKK.GE. 123. AND.KK.LE. 146) BETA (KK) =3. 14136 
IFCKK.GE. 147. AND.KK.LE. 158) BETA(KK)=3 . 29843 
IFCKK.GE. 159 .AND.KK.LE. 202) BETA (KK) =3 .40314 
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IF (KK . GE . 203 . AND . KK . LE . 2 14) 
IFCKK.GE. 215, AND.KK.LE. 230) 
IFCKK.GE.231.AND.KK.LE.246) 
IF (KK . GE , 247 . AND . KK . LE . 25 8 ) 
IFCKK.GE .259 . AND.KK.LE . 326) 
IFCKK.GE.327.AND.KK.LE.350) 
IFCKK.GE.351.AMD.KK.LE.430) 
IFCKK.GE .431 .AND. KK. LE .438) 
IFCKK.GE .439 . AND.KK.LE .446) 
IFCKK.GE. 447 .AND.KK.LE.454) 
IFCKK.GE. 455. AND.KK.LE. 462) 
IFCKK.GE. 463. AND.KK.LE. 470) 
IFCKK.GE. 471. AND.KK.LE. 478) 
IFCKK.GE. 479. AND.KK.LE. 486) 
IFCKK.GE. 487. AND.KK.LE. 494) 
IFCKK.GE. 495. AND.KK.LE. 502) 
IFCKK.GE. 503. AND.KK.LE. 510) 
IFCKK.GE. 511. AND.KK.LE. 518) 
IFCKK.GE. 519. AND.KK.LE. 526) 
IFCKK.GE. 527. AND.KK.LE. 534) 
IFCKK.GE. 535. AND.KK.LE. 542) 
IFCKK.GE. 543. AND.KK.LE. 550) 
IFCKK.GE. 551. AND.KK.LE. 558) 
IFCKK.GE. 559. AND.KK.LE. 638) 
IFCKK.GE. 639. AND.KK.LE. 662) 
IFCKK.GE. 663. AND.KK.LE. 708) 
IFCKK.GE. 709. AND.KK.LE. 752) 


BETACKK)=3.29843 
BETA (KK)=3. 14136 
BETA(KK)=0.0 
BETA CKK)=-. 15707 
BETACKK)=-.2618 
BETA(KK)=“. 15707 
BETA CKK) =0.0 
BETACKK)=. 36394 
BETA (KK)=. 17452 
BETA(KK)=. 38394 
BETACKK)=. 17452 
BETA(KK)=1 . 2216 
BETACKK)=. 69808 
BETA(KK)=1 . 2216 
BETA(KK)=. 69808 
BETA CKK)=2. 4433 
BETA(KK)=1. 91972 
BETA CKK )=2. 4433 
BETACKK)=1 .91972 
BETACKK)=2.9668 
BETA(KK)=2. 7574 
BETA(KK)=2.9668 
BETACKK)=2.7574 
BETACKK)=3. 14136 
BETACKK)=3. 29843 
BETACKK)=3. 40314 
BETACKK)=0.0 


READCINj 1) NEL,NS > XC1) > YC1) J XC2) J YC2) > XC3) 3 Y(3) 

C VRITEC6 ,121) NEL J NS J XC1) ) YC1) J XC2) J YC2),XC3),YC3),BETACKK) 
121 F0RMATC4I5 , 7F10.4) 


DIMENSION IN FEET... 

X1=XC1) 

X2=XC2) 

X3=XC3) 

Y1=YC1) 

Y2=YC2) 

Y3=YC3) 

C 

X C 1 5= X(l) /12.0D0 
X ( 2 )= XC2) /12.0D0 
X ( 3 )= XC3) /12.0D0 

Y C 1 )= YC1) /12.0DO 

Y C 2 )= YC2) /12.0D0 

Y ( 3 )= YC3) / 12. ODD 
C 

T1=T(1) 

T2=T(2) 

T3=T(3) 

1107 FORMAT C1H1,///20X 

* , ' TCI) TC2) TC3) 

* ISIDEC1) ISIDEC2) ') 

1108 FORMAT C' ' ,20X 
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t _ 

* / ') 

G WRITE (10, 11)T, ISIDEH(KK, 1) , ISIDEH(KK, 2) 

C 

C 

c 

T0LD(NS(1))=T(1) 

T0LD(NS(2))=T(2) 

TOLD(NS(3))=T(3) 

C H : CONVECTION COEFFICIENT 

C WRITE (6 ,5152) KK,NS(1) ,NS (2) ,NS(3) ,T0LD(NS(1) ) ,T0LD(NS(2) ) , 

C 1 T0LDCNSC3)), H(KK),(ISIDE(JJ),JJ=1,2), 

C 1 XI ,Y1,X2,Y2,X3,Y3,NGT(KK, 1) ,N0T(KK,6) 

5152 FORMAT C' ’ ,415 ,3F10. 3 ,F5 .0 , 212 , 6F8 .4 , 15 , AS) 

C 

T(1)-T1 

T(2)=T2 

T(3)=T3 

C 

C ECHO THE INPUT 

C WRITE (10, 23) NEL,NS ,X(1),Y(1),X(2),Y(2),X(3),Y(3) 

C * , HH,RHO,CP, ISIDE ,T 

23 F0RMAT(1H0, 15 , 1X,3I5 , 2X,6(1X,F8.4) ,2X, 2F5 . 0,F6 . 3,2X,2I3, 

* 2X,3F6.Q) 

C CALCULATION OF THE CONDUCTION MATRIX 
C 

C KPT (NE ,3) : ARRAY STORES THE NODE NUMBERS FOR EACH ELEMENT 
NPT(KK, 1) = NS(1) 

NPT(KK,2) = NS (2) 

NPT(KK,3) = NS (3) 

RCNS(l)) = X(l) 

R(NS(2) ) = X(2) 

R(NS(3)) = X(3) 

Z(NSCD) =Y(1) 

Z (NS ( 2 ) ) — Y (2) 

Z(NS (3) ) =Y(3) 

AREAE(KK)=(X(2)*Y(3)-X(2)*Y(1)+X(3)*Y(1)-X(3)*Y(2) 

1 + X(1)*Y (2) -X(1)*Y(3) )/2, 0 


C C! -C 

C - 

C WRITE (6 ,266)KK,X(1) ,Y(1) ,X(2) ,Y(2) ,X(3) ,Y(3) ,AREAE(KK) 

266 FORMAT (1H0 , 15 , 6(1X,F8 .4, ) ,2X,D14. 7) 

C WRITE(6,319)KK,(T(I) ,1=1,3) 

319 FORMAT( t ELEMENT NO. ' ,I8,3F12.4) 

7 CONTINUE 

C END OF THE DO LOOP Q\ 

C WRITE (6, 264) (I, TOLD (I) , 1=1 ,NP) 

264 FORMAT (1H0, 15 ,E14.5 ,2X,I5 ,E14.5 ,2X, 15 ,E14.5 , 2X, 15 ,E14.5 


* , 2X,I5,E14.5,2X,I5,E14.5) 

C 

C ADD EXTERNALLY APPL. CONC, NODAL FORCES VECTOR TO EF VECTOR 
ID 1=0 
INK=0 

202 READ (IN, 203) NODNUM,CFV 
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203 FORMAT ( 6 I 3 , 2X , 6F10 . 5 ) 

ID=0 

DO 204 L-l , 6 

IF C NODNUM(L).LE.O) GOTO 20S 

ID=ID-f-l 

I NQDHU=NODNUM ( L ) 

204 EF ( INODNU) =CFV (L)+EF ( INODNU) 

GOTO 206 

205 INK=1 
IF(ID.EQ.O) GOTO 216 

206 IF(ID1 .EQ. 1)G0T0 222 
WRITE (10, 208) 

208 FORMAT (1H1, ' NODE NUMBER ',10X,’ GENERALIZED NODAL FORCES ') 
222 WRITE(IOj207) (NODNUM(L) ,CFV(L) ,1=1,10) 

207 FORMAT (1H0, 6(13 ,E 14. 5 ,2X) ) 

IF(INK.EQ.1)GOTO 216 
ID1=1 

GOTO 202 
216 CONTINUE 

C 

JJJKNT=0 

C 

DO 3000 ITER=1, 100 

CALL FORMTK ( NE , NP , NBW , TIME , ITER ) 

CALL FORMTC (NE,NP,NBW) 

CALL FORMTF (NE,NF,NBW) 

CALL ASSMBL (NP, NBW, DDT, ITER) 

C WRITE (6 , 9999) ( (ZZ(I , J) , J=1,NBW) ,1=1, NP) 

9999 F0RMAT(2X,7G11.5) 

C WRITE (6, 9988 ) (F (1,1), 1=1, NP) 

9988 FORMAT (2X, Gil. 5) 

CALL SOLVE (NP, NBW) 

JJJKNT=JJJKNT+l 
IFfJJJKNT.GE. 1)G0TQ 979 
GOTO 980 

979 CONTINUE 

TIME ( ITER ) =DDT*ITER 
DO 2001 I = 1,483 
FI J( I, ITER) = F(I) 

2001 CONTINUE 

C WRITE ( 10 , 666)TIME , ITER , TDT 

C66 FORMAT(/ / , 1H1 , 1 TEMPERATURE DISTRIBUTION AT THE TIME =' , 

C * E15.6, ' HOUR ' ,/' ITERATION NO. =’,110, , 

C * 1 *&* TDHEAT . 2 DDT= r ,F10.6,/) 

C IF(ITER.EQ. 22) 

C 1 WRITE (7, 777) (I,F(I) ,I=1,NP) 

C77 FORMAT C 1H0 ,15 ,F14. 5 , 2X, 15 ,F14.5 ,2X, 15 ,F14.5 , 2X, 15 ,F14. 5 , 

C *2X, 15 ,F14.5 ,2X, 15 ,F14. 5) 

C77 FORMAT (I4,F10.4) 

980 CONTINUE 
MAX(1)=0 . 0 

srr=o . o 

DO 1515 1=1 ,NP 


A 'll 



CF 

1211 


1515 

C 

CF 

CF 

1222 


2000 

3000 

1200 

666 

C 

777 

3011 

C 

C 

C200 


** « 
c.. 
c 


c 


c 

c 

c 

c. 


MAX(I+1)=DABS CF( X) -TOLD ( I ) ) 

WRITE (6, 1211 )MAX(I+1) 

FORMAT (2X,F 15 . 6) 

IF (MAX (1+1) .LT.MAX(I) )GOTO 1515 
IF (MAX (1+1 ) .GT.MAX(I) ) FAXM=MAX (1+1) 

IF(FAXM . GT . STT) GTST-FAXM 

STT=GTST 

CONTINUE 

WRITE (6 , 1222) GTST 

GTST IS THE STOPPING CRITERION. ITERATION STOPS WHEN DIFF. BETW 
T(NEW) AND T(OLD) IS LESS THAN 0.3 DEGREE 
FORMATC/ ,2X,F10.5) 

IF (GTST. LT. 0.3) GOTO 1200 
DO 2000 1=1 ,NP 
TOLD(I)=F(I) 

CONTINUE 

CONTINUE 

DO 3011 IJP = 1,483 
WRITE (7 ,666) IJP 

FORMAT(//, t TEMPRATUE DISTRIBUTION AT POINT ' ,15,/) 

WRITE (7,777) (I ,TIME(I) ,FIJ (IJP, I) , I =1,65) 

FORMAT (15 ,2F10 . 4) 

CONTINUE 
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STOP 

END 

SUBROUTINE FORMTK (NE,NP,NBW, TIME, ITER) 






FORM BANDED THERMAL CONDUCTANCE MATRIX 


IMPLICIT REAL*4(A-H,0-Z) 

REAL'-' 8 KZETA , KETA , RHOC P , BETA 

REAL*8 BI , BJ, BK, Cl , CJ, CK, COF , CONI ,C0N2 , C0N3 

REAL* 8 REAR, CON , LG, RCTL ,ELKT( 3 ,3) ,THSTF (483 , 18 ) 

REAL*8 AGO, All , A22,BQ0 ,B11 ,B22,CQ0,THCAP(4B3, 18) ,TOLD(483) 
REAL*8 Z(483) ,R(483) ,AREAE (752) ,X(3) ,Y(3) , DDT, TIME, ATIM,P 
COMON/CONVEN/H(752) ,ISIDEH(752,2) 

COMMON / SOLID / Z , R , AREAE , THSTF 
1 ,THCAP,NPT(752, 3) ,T0LD 

REAL*8 ZZ(483, 18) ,F,B,FF,EF 
COMMON / SOLVEQ / F(483) ,EF(752) ,ZZ 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RH0CP(752) ,BETA(752) 


INITIALIZE GLOBAL THERMAL STIFFNESS MATRIX 

IN=5 

10=6 


DO 120 J=1,NBW 
DO 110 1=1 ,NP 
THSTF (I, J) = O.D0 
110 CONTINUE 


A Zfc 


o o 


120 CONTINUE 
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c 

c 

130 NBAND = 0 

DO 899 J=1,NE 

C WRITE (6, 788) KZETA (J ) ,KETA(J) ,BETA(J) 

788 FORMAT (3X,3F20. 9) 

899 CONTINUE 
C 

C WRITE (6,7776) 

7776 FORMAT ( * CONDUCTION MATRIX ELEMENTS ARE ',/) 

DO 230 N-1,NE 

Coooooo <><><><><><><><><><><> 

C 

II = NPTCN.l) 

JJ - NPTCN.2) 

KK = NFT(Nj 3) 

BI - Z(JJ)-Z(KK) 

BJ = Z(KK) -Z(II) 

BK = Z(II)-Z(JJ) 

Cl = R(KK)-R(JJ) 

CJ - R(II)-R(KK) 

CK “ R(JJ)-R(II) 

X(l) ~R(II) 

X(2) -R(JJ) 

X(3) =R(KK) 

Y(l) =Z(II) 

Y(2) =Z(JJ) 

Y{3) =Z (KK) 

C WRITE(6 j555)X(l) ,X(2) ,X(3) ,Y(1) , Y(2) ,Y(3) 

555 FORMAT ( 2X , 6F 12 . 6 ) 

ALPHA=BETA(N) 

C0N1=(KZETA(N)* (COS (ALPHA) )**2+KETA(N)*(SlN (ALPHA) )**2) 
1/(4.*AREAE(N)) 

CON 2= (KZETA (N) -KETA(N))*SIN(2.*ALPKA) / (4.*AREAE(N)) 

CON3= (KZETA (N) * ( S IN ( ALPHA ) )**2+KETA(N)*(C0S (ALPHA) )**2) 
1/(4.*AREAE (N)) 

WRITE (6,7999) CONI , C0N2 , C0N3 , AREAE (N) 

7999 FORMAT(2X,3F10.3,2X,F10. 7) 

ELKT(1, 1) = ( B 1**2 )*C0N1+ ( C 1**2 ) *CON3 
ELKT (1,2) = (BI*B J)*C0N1+ (CI*CJ)*C0N3 
ELKT(1,3) = (BI*BK)*C0N1+(CI*CK)*C0N3 
ELKT (2, 2) = ( B J**2 )*C0N 1+ ( C J* *2 ) *CON3 
ELKT(2,3) = ( B J*BK) *GONl+ ( CJ--CK ) * COM3 
ELKT(3,3) = (BK**2) *C0Nl-f-(CK ,v *2 ) *CON3 
DO 10 IQ=1,2 

C WRITE(6,7805)N, ( ISIDEH(N, Jl) , Jl=l,2) 

7805 FORMAT ( 1 SIDE EXPOSED ’,3110) 

IF ( ISIDEH(N, IQ) . LE . 0 ) GO TO 240 
JQ = ISIDEH (N, IQ) 

C WRITE (10 ,12) JQ,N 

12 FORMAT (1H0 , 1 INSIDE* ,13, 3X, * OF ELEMENT* ,15) 

KQ = JQ +1 
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IF ( JQ.EQ . 3) KQ=1 

LG s DSQRT(CXCKQ)-XCJQ))**2+CYCKQ)-Y(JQ))— 2) 
HL =H(N)*L6/6. 

IF(JQ.EQ.l) GO TO 20 
IFCJQ.EQ.2) GO TO 30 
IF(JQ.EQ.3) GO TO 40 
GO TO 10 

20 ELKT(1 j 1)= ELKT(1,1)+HL * 2. 

ELKT(2 ,2)= ELKT(2,2)+HL * 2. 

ELKT(1,2)= ELKT(1,2)+HL 
GO TO 10 
C 

30 ELKT(2,2)= ELKT(2, 2)-rHL*2 . 

ELKT(3 ,3)= ELKT(3,3)+HL'‘2. 

ELKT(2,3)= ELKT (2,3) +HL 
GO TO 10 
C 

40 ELKTCl s 1)= ELKT(1 , 1)+HL * 2. 

ELKT(3,3)= ELKT (3,3) +HL * 2. 

ELKT(1 ,3)= ELKT(1 3 3)+HL 
10 CONTINUE 


WRITE CIO, 2000) N, II, JJ,KK,XC1) ,Y(1) ,X(2) ,Y(2) ,X(3) * Y(3) ,CONDTY(N) , 
1AREAECN) I RBAR,CQN,ELKTC1 J 1),ELKTC1,2),ELKTC1,3),ELKTC2,2), 
2ELKTC2,3),ELKTC3,3) 

2000 FQRMATC 1HG > 4 IS , 8G11 . 4 a /21X , 8G11 . 4) 

— - - -- — - --C! 

STORE IN GLOBAL MATRIX CS 

— — — C! 

240 GONTINUE 

WRITE C IO,266)N,X(l) ,Y(1) ,X(2) ,Y(2) ,X(3) ,Y(3) , AREAE(N) ,NBAND 
266 FORMAT C1H0 , IS , 6 C IX , F8 - 4) , 2X, D14 . 7 , 2X , 15 ) 

WRITE (6 j 277 )N , C( ELKTCl, J) ,J=1,3) ,1=1,3) 

277 FORMAT C * ',15/ 3(3F25.6/)) 

DO 220 LL=1 , 3 
DO 210 MM-1,3 
IF(MM.LT.LL) GOTO 210 
I=NPT(N,LL) 

J=IABS CNPT CN , MM)-I)+1 

IFCNPTCNjMM) .LT.NPT(N,LL)) I=NPT(N,MM) 

THSTF C I , J) =THSTF C I , J)+ELKT CLL , MM) 

WRITE (6, 7543)1, J, II ,JH,NPT(NjJ) ,NPT(N, I) 

7543 FORMAT C ' ' ,613, '******•;) 

IF (J.GT.NBAND) NBAND=J 
210 CONTINUE 
220 CONTINUE 
230 CONTINUE 

WRITE C6 , 5 05 0 ) C (THSTF ( I , J) , J=1 , NBW) , 1=15 , 20 ) 

5050 FORMAT(7G11.5) 

WRITE (6,200) NBAND 
200 FORMAT C1H0 , 1 NBAND= ' ,15) 


END OF PROGRAM 
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c 

RETURN 

END 

C 

C FORM BANDED CAPACITANCE MATRIX FOR A DEFORMING F.E. MESH 

C 

SUBROUTINE FORMTC (NE.NP.NBW) 

C 

REAL*8 ZC483) ,R(483) ,AREAE(752) ,X(3) ,¥(3) ,DDT 
REAL* 8 CON t ELCAP(3,3) ,THCAP(483 , 18) ,T0LD(483) ,THSTF(483 , 18) 
1 ,THSTF1(747,17) 

COMMON/CONVEN/H ( 75 2 ) , S IDEH (752,2) 

C 

COMMON / SOLID / Z , R , AREAE , THSTF 
1 , THCAP , NPT ( 75 2 , 3 ) , TOLD 
READ'S ZZ (483 , 18 ) , F , EF , B , COF , FF , KZETA , KETA 
COMMON / SOLVEQ / F(483) ,EF(752) , ZZ 
DIMENSION FF( 15 , 1) 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RHQCP(752) , BETA(752) 

C 

C 

REAL * 8 R1,R2,R3,R4,R5,R6 
C 

c INITIALIZE GLOBAL THERMAL CAPACITANCE MATRIX 

C 

10=6 

J9=0 

DO 120 J=1,NBW 
DO 110 1=1 ,NP 
THCAP(I,J) = O.DO 
110 CONTINUE 
120 CONTINUE 
C 
C 


C WRITE (6,8880) 

8880 FORMAT (‘ CAPACITANCE MATRIX ELEMENTS ARE 1 ,/) 

200 DO 230 N=1,NS 
C 

G FOR EACH ELEMENT, FOR THERMAL CAPACITCANCE MATRIX 


G 

II = NPT(NYl) 

JJ = NPT(N, 2) 

KR = NPT (N } 3 ) 

CON=RHOCP (N) * AREAE (N) / 12 . 

ELCAF ( 1 , 1 ) =2 . *COM 
ELCAP(1,2)=1. *CON 
ELCAP ( 1 , 3 )=1 , *CON 
ELC AP ( 2 , 2 ) =2 . *CON 
ELC AP ( 2 , 3 ) =1 „ *CON 
ELC AP ( 3 , 3 ) =2 - *CON 
C 
C 

C WRITE(IQ,2000)N, II, JJ,KK,R(II) ,Z(II) ,R(JJ) , Z(JJ) ,R(KK) , Z(KK) , 



C 1 RHOCP(N) ,AREAE(N) ,CON,ELCAP(l, 1) ,ELCAP(1,2) ,ELCAP(1,3) , 
G 2 ELCAP(2,2),ELCAP(2,3),ELCAPC3,3) 

2000 FORMAT ( XHO , 415 , 9G1 1 . 4/2 IX , 6G11 . 4) 

C 

C STORE IN GLOBAL MATRIX 


DO 220 1=1,3 
DO 210 J=l,3 
IF(J . LT. I) GOTO 210 
II=NPT(N,I) 

JJ=IABS(NPT(N\J)-II)+1 
IF (NPT(N,J).LT.NPT(N,I)) II=NFT(N,J) 
THCAP ( II, JJ) « THCAP(II,JJ)+ELCAP(I,J) 


210 CONTINUE 

220 CONTINUE 

C 

C WRITE (6 , 2 7 7 )N , ( C ELCAP(I , J) , J=1 ,3) , 1=1, 3) 

277 FORMATf ' , 15/ 3(3F25 . 6/) ) 

230 CONTINUE 
C 

C WRITE (10,4000) CTHCAP ( I , J) , J=1 , NBW) 

4000 • FORMAT (1H0 , 1 THCAP ’ , 10GU . 4/6X , 10G1 1 . 4/6X , 10G11 . 4/6X , 10G11 . 4) 
C WRITE (10, 1000) 

1000 FORMAT C1H0 , ' LEAVING FORMTC IN TDHEAT ..IB ! ) 

264 FORMAT (1H0 , 15 ,E14. 5 ,2X, 15 ,E14. 5 ,2X,I5 ,E14. 5 , 2X, 15 ,E14.5 

* ,2X, 15 ,E14.5 , 2X, 15 ,E14.5) 

C 

C END OF PROGRAM 


C 

RETURN 

END 

SUBROUTINE FORMTF (NE , NP , NBW) 

C 

C FORM THERMAL FORCE VECTOR 

C 
C 
C 

REAL*8 ELF (3) , LG , RCTL , COF , THCAP (4B3 , IB ) , THSTF (483,18) 

REAL*S 2(483) ,R(483) S AREAE (752) ,X(3),Y(3) , DDT, TOLD (48 3) 

COMMON/ CONVEN/ H(752),ISIDEH(752,2) 

C 

COMMON / SOLID / Z ,R, ARE AE, THSTF 
1 , THCAP, NPTC752, 3), TOLD 

REAL* 8 Z Z ( 48 3 , 1 8 ) , F , EF , B , E LF 1 , ELF2 , ELF3 , FF , HLL 

COMMON / SOLVEQ / F(483) ,EF(752) , ZZ 
REAL *8 KZETA , KETA , RHOCP , BETA 

COMMON /CONT/ KZETA(752) ,KETA(752) ,RH0CP(752) , BETA(752) 

REAL *8 QR(752) ,HL(752) 

DATA W / . 6666666666D0/ 

10=6 

TINF =10. DO 
C 
C 

C FOR EACH ELEMENT , FORM THERMAL FORCE VECTOR 

C 

A 32 


c 

9999 

110 

1070 


C 

999 

20 

30 
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WRITE C6 S 9999) 

FORMAT C FORCE VECTORS ARE ’,/) 

DO 110 1=1 ,NE 

EF(I) = 0 .DO 

QRCIJ^O.O 

CONTINUE 

DO 1070 J=l,45 

N=2*J 

QR(N)=326.0 

CONTINUE 

DO 230 N=1 ,NE 

ELF(1)=0.D0 

ELF(2)=0.D0 

ELF(3)=Q.D0 

II=NPT(N, 1) 

JJ=NFT(N,2) 

KK=NPT(N,3) 

Cl = RCKK)-R(JJ) 

CJ = R(II)-R(KK) 

CK = RCJJ)-R(II) 

BI=Z(JJ)-ZCKK) 

BJ=ZCKK)-Z(II) 

BK=Z(II)-Z(JJ) 

X(1)=R(II) 

X(2)=RCJJ) 

X(3)=RCKK) 

YC1)«ZCII) 

YC2)=ZCJJ) 

YC3)=Z(KK) 

ALPHA=BETA(N) 

C0N2=(KZETA (N) -KETA(N) )*SIN (2 .* ALPHA) / (4 . *AREAE (N) ) 

ELF 1=C0N2* CB I*CI*T0LD (II) +BI*C J*TOLD C JJ)+B I*CK*T0LD (KK) ) 

ELF2=C0N2* (BJ*CI*TOLD ( I I)+BJ*CJ*TOLD (JJ)+BJ*CK*TQLD (KK) ) 

ELF3=C0N2*(BK*CI*T0LD (II) +BK*C J*TOLD ( J J) +BK*CK*TOLD (KK) ) 

ELF(1)=~W*ELF1 

ELF (2)=“W*ELF2 

ELF (3 )=“W*ELF3 

DO 10 10=1,2 

IF(ISIDEH(N,IQ).LE.O)GO TO 240 
JQ=ISIDEH (N, IQ) 

KQ=JQ+1 

IFCJQ.EQ.3) KQ=1 

LG=DSQRT( (X(KQ) -X( JQ) )**2+ (YCKQ) -Y( JQ) )**2) 

WRITE C6,999)LG 
F0RMAT(2X,F10.5) 

HL (N ) =QR (N) +H (N ) *TINF 
IF ( JQ-EQ. 1)G0 TO 20 
IF C JQ-EQ. 2)G0 TO 30 
IF (JQ .EQ.3)G0 TO 40 
GO TO 240 

ELF (1) =HL CN ) *LG* . 5+E LF ( 1 ) 

ELF (2) =HL (N) ‘"LG* . 5+ELF ( 2 ) 

GO TO 10 

ELF (2) =HL(N)*LG*. 5+ELF C2) 


on o o ooh-^oonoo 
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ELF (3) =HL(N)*LG*.5+ELF(3) 

GO TO 10 

40 ELF(l) =HL CN ) *LG* f 5+E LF ( 1 ) 

ELF(3) =HL(N)*LG*.5+ELF(3) 

10 CONTINUE 

C 

240 EF(II)= EF(II)+ELF(1) 

EF(JJ)=EF(JJ) +ELF(2) 

EF(KK)=EF(KK) +ELF(3) 

707 FORMATC 1 ,2I8,5F15.5) 

IFCEF(II) .EQ.O) GO TO 50 
IFCEF(JJ).EQ.O) GO TO 50 
IFCEF(KK).EQ.O) GO TO SO 
GOTO 230 

WRITE (6, 1999)N,II , JJ,KK,EF(II) ,EF(JJ) ,EF(KK) 
FORMATC 1H0 ,416 , 2X, 3F10 .4) 


230 

CONTINUE 

WRITE (10, 2000) (EF (I), 1=1 ,NP) 


2000 

FORMAT (2X.F10. 4) 
WRITE (10, 1000) 


1000 

FORMAT ( 1H0 , ' LEAVING FORMTF 

RETURN 

END 

IN TDHEAT . . .IB') 


SUBROUTINE ASSMBL (NP,NBW, DDT, ITER) 

ASSEMBLE MATRICES FOR RECURRENCE FORMULAS 

REAL*8 Z(483) ,RC483) ,AREAE(752) ,THSTF(483 , 18) ,THCAPC483 , 18) 
1 ,THSTF1(747 , 17) 

REAL*8 ZZ(483,lS),F,EF,B,ZO(747,17),FF(15,l) s TOLD(483) 

COMMON/ CON VEN/H ( 75 2 ) , IS IDEH (752,2) 

COMMON / SOLID / Z , R , AREAE , THSTF 
1 ,THCAP,NPT(752 , 3) ,TOLD 
COMMON / SOLVEQ / F(483) ,EF(752) ,ZZ 

COMMON /CONT/ KZETAC752) ,KETA(752) ,RH0CP(752) ,BETA(752) 

WEIGHTING FACTORS FOR GALERKIN ANALYSIS 


DATA W1,W2,W3,W4/ 0.333333333333333D0 , 

1 0.500000000000000DO , 

2 “0, 1666666666666667D0 , 

3 0.500000000000000DO / 

C 

10=6 

DT = DDT 

C WRITE (10,1000) W1,W2,W3,W4,DT 

1000 FORMAT (1H0, 'ASSEMBLE IN TDHEAT 1 1 // ,5E15 . 8/ T ZZ 

1 COL 1 TO 6 '/) 

C 

6567 FORMATC ',I5,6F12.6) 

3000 FORMAT (1H0 f ' THSTF ' , 10G11 . 4/6X, 10G1 1 . 4/6X , 10G1 1 . 4/6X, 10G1 1 . 4/ ) 
4000 FORMAT (1H0 , ' THCAP ' , 10G11 . 4/6X , 10G11 . 4/6X , 10G1 1 . 4/ 6X , 10G1 1 . 4/ ) 
C 

DO 150 1=1, NP 
C 


A 34 
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COEFFICIENT MATRIX 

DO 110 J-1,NBW 

ZZ(I,J) = W1*(THSTF(I,J))4W2*CTHCAP(I,J))/DT 
110 CONTINUE 

WRITE (10 , 2000) (22 ( I , J) , J-l , NBW) 

2000 FORMAT (' ' ,8E15.6) 

THERMAL LOAD VECTOR 

F C I ) - EFCD/2.D0 

IF (I.EQ.l) GO TO 130 
JST = I+l-NBW 
IF (JST.LT.l) JST=1 
JEND = 1-1 

WRITE (6,1000) I, JST, JEND 

DO 120 J=JST, JEND 
II = I+l-J 

- 01 — - 

F ( I ) =F(I)+ (W3* (THSTF ( J , 1 1 ) ) +W4* (THCAP ( J , 1 1 ) ) 

1 /DT)* (TOLD(J)) 

WRITE (6,1000) I , J, II ,THSTF(J, II) , THCAP (J, II) ,TOLD(J) ,F(I) 
120 CONTINUE 

130 JST = I 

JEND = I-l+NBW 
IF (JEND . GT . NP ) JEND=NP 

WRITE (6,1000) JST, JEND 
DO 140 J= JST, JEND 
JJ = J-I+l 

F ( I ) =F( I)+(W3* (THSTF ( I , JJ) )+W4* (THCAP ( I , JJ) ) 

1 /DT)* (TOLD(J)) 

WRITE (6,1000) I,J,JJ,THRTF(I,JJ),THCAP(I,JJ),TOLD(J),F(I) 
140 CONTINUE 

WRITE (6,1000) I , JST, JEND ,F(I) 

WRITE (6,1001) I, (ZZ(I,J) , J— 1,NBW),F(I) ,TOLD(I) 

150 CONTINUE 

1001 FORMAT (1H ,I3,6G11.5) 

END OF PROGRAM 

RETURN 

END 

SUBROUTINE SOLVE (NP, NBW) 

C 

C SOLUTION OF EQUATIONS BY GAUSS ELIMINATION 

C NP=NO. OF NODAL POINTS NBW=B ANDW IDTH 

C =NO. OF UNKNOWNS 

C =NO. OF EQUATIONS 

C 

REAL*8 ZZ (483, 18), F,EF,B,P, EXTRA 
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C 

1002 

100 

c 


1001 

c 

c 


c 


110 


c 

c 

c 

c 

120 


c 

c 


1000 

130 

c 

c 

c 

200 

210 


220 


C 

C .. 
C 

300 


COMMON/SOLVEQ/F(483) ,EF(752) ,ZZ 
10=6 

WRITE (6, 1002) 

FORMAT (* IN SOLVE ZZCI.J) COL 1 TO 4 '/) 

DO 130 1=1, NP 

WRITE (10, 1001)1 , (ZZ(I , J) , J=l,4) ,F(I) 

FORMAT (* ' ,I5,5E15.5) 

REDUCE THE STIFFNESS MATRIX 

DO 120 J=2,NBW 
II=I+J-1 

IF(ZZ(I,J).EQ.O.O) GO TO 120 
P= ZZ(I,J)/ZZ(I,1) 

JJ=0 

DO 110 K=J , NBW 
JJ=JJ+1 

IF(ZZ(I,K).NE.O.ODO) ZZ (II , JJ)=2Z (II , JJ) -P*ZZ(I ,K) 

IF(DABS(ZZ(II,JJ)).LT. l.E-20)ZZ(II, JJ)=0.0D0 
CONTINUE 
ZZCI,J) «P 

REDUCE LOAD VECTOR 
F(II) = F(II)-P*F(I) 

CONTINUE 

F(I) = F(I)/ZZ(I , 1) 

WRITECIO, 1000)1, ZZ(I,J),J=1, NBW), F(I) 

FORMAT ( 1H0 , ’ SOLVE * , 1 3/4 ( 5 E IS . 8 / ) ) 

CONTINUE 

BACK SUBSTITUTION 


N=NP 
N=N- 1 

IF(N.EQ.O) GO TO 300 
II=N 

DO 220 J=2,NBW 
11 = 11+1 

IF ( ZZ(N, J) .NE.O.O) F(N)=F(N)"ZZ(N, J)*F(II) 

CONTINUE 

GO TO 210 

.SOLUTION VECTOR STORED IN F ARRAY 

RETURN 

END 
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c 

C THE GRID PROGRAM IS MODIFIED TO GENERATE RECTANGULAR ELEMENTS, 

C THE CHANGES ARE AS FOLLOW: 

C 

CH THE DIMENSION OF N00C2500,6) IS CHANGED TO N00(2500,7) 

CH THE DIMENSION OF NR(3) IS CHANGED TO NR (4) 

CH THE DIMENSION OF IS CHANGED TO LB (3) TO LB (6) 

C 

C DIVISION INTO RECTANGULAR ELEMENTS 
C 

K-l 

DO 17 1=1 ,NROWS 
DO 17 J=1 ,NCOL 
XE(K)=XC(I , J) 

YE(K)=YC(I,J) 

NE(K)=NN(I,J) 

17 K=K+1 
L=NROWS-l 
DO 21 1=1, L 
DO 21 J=2 , NCOL 

GH DIAS1=SQRTC(XCCI,J)-XCCX+1,J-1D)**2+(YCCI.J)“YC(I+1,J-1))**2) 

CH DIAG2=SQRT ( (XC C 1+1 , J) -XC ( I , J-l) )**2+ ( YC ( 1+1 , J) - YC ( I , J - 1 ) )**2 ) 

NR ( 1 ) -NCOL* 1+ J - l 
NR ( 2 ) =NCOL* 1+ J 
NR ( 3 ) =NCOL* ( I - 1 ) +J 
NR(4)=NC0L*(I-l)-hJ-l 
CH DO 21 IJ=1,2 

NEL=NEL+1 

CH IFC CDIAG1/DIAG2) . GT. 1. 02) Stt TO 18 

J1=NR(1) 

J2=NR(2) 

J3=NR(3) 

J4 = NR (4) 

CH GO TO 19 

CH IS J1=NR(IJ) 

CH J2=NR(IJ+1) 

CH J3=NR(4I 

CH LB (3) ,LB(5) AND LBC6) ARE ADDED 
LB C 1)=IABS CNE ( J1 ) -NE ( J2 ) )+l 
LBC2)=IABSCNE(J1)-NE(J3})+1 
LB C 3 ) =1 AB S CNE ( J1 ) -NE ( J4 ) ) +1 
LB(4)=IABS(NE(J2)-NE(J3) )+l 
LB (5 )=IABS (NE ( J2) -NE ( J4) )+l 
LB C6 )=IABS CNE ( J3) -NE ( J4) )+l 
DO 20 IK=I,6 

IFCLB(IK) .LE.NBW) GO TO 20 
NBW=LBCIK) 

NELBW=NEL 
20 CONTINUE 

IFC KK . NE . 46 ) GO TO 688 
C WRITEC6,666) J3 , NECJ3) 

666 FORMAT C///// ' J3 = ’,19,' NE(J3) = ' , 19 ) 

888 CONTINUE 

WRITE CIO, 113) KK,NEL,NE(J1),NE(J2),NE(J3) ,NE CJ4) ,XE CJ1) ,YE C Jl) , 
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1 XECJ2) , YE(J2) ,XECJ3) ,YE(J3) ,XE(J4) , YE(J4) 

IFCIPCH.EQ.O) GO TO 21 

WRITE C IP , 1 14 ) NEL , NE C J1 ) , NE ( J2 ) , NE ( J3 ) , NE C J4) , XE ( J 1 ) , YE ( J1 ) , XE ( J2 ) 
1,YE(J2) ,XE(J3) ,YECJ3),XECJ4) ,YE(J4) 

N00 CNELjfi) IS ADDED 
NOOCNEL, 1)=KK 
NOO CNEL, 2)=NEL 
NOOCNEL, 3)=NE(J1) 

N00CNEL,4)=NE(J2) 

N00CNEL,5)=NECJ3) 

NOOCNEL, 6 )=NE(J4) 

NOOCNEL, 7 )=IRRCKK) 

21 CONTINUE 

22 CONTINUE 
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CH THE DIMENSION OF N00, X, I, OR Y IS INCREASED BY 1 
C 

COMMON/AAA/ N00(2500,7) ,NTT(80) 

COMMON/BBB/ X(4,2000) ,Y(4,2000) ,NBO(1000) S TITLEC10) 

COMMON/ CCC/ NODES , LMENTS ,JT (12000), MEM JT C 24000 } , JMEM (3000) , 
$JNT(3000) , IDIFF,MINMAX 
DATA IN/11/ , 10/6/ , IF/12/ 

C 

READ (IN, 100)TITLE 
WRITE (IP, 100)TITLE 
100 FORMAT (10A4) 

122 FORMAT ( 13 jSFIO .5) 

J =1 
JJ=0 

CH I IS CHANGED FROM 3 TO 4 

150 READ(IN, 17) ( JT(3000*(I-l)'hJ) , 1—1,4) , (X(I J) , Y(I , J) ,1=1 ,4) 

CH JT(9000+J)=0 

C WRITE(IO,17)CJT(3DOO*(I-1)+J),I=1,4),(X(I,J),YCI,J),I=1,4) 

NODES=MAXO(JT(J) , JT(3000+J) , JT(6000+J) , JT(9000+J) , JJ) 
JJ=N0DES 

IF(JT(3000+J).EQ.O)GO TO 152 
J=Jt1 
GO TO 150 
152 LMENTS=J-1 

NODES=JJ 

12 FORMAT (/ / , 5X , 15HNUMBER OF NODES, 14, 15X, 

1 T NUMBER OF ELEMENTS ',14,//) 

WRITE CIO, 105 )TITLE 
105 FORMATC//, 1X,20A4,//) 

C WRITE ( 10, 12)N0DES, LMENTS 

C WRITE (10, 13) 

C3 FORMAT (//,3X,18HNEL JT1 JT2 JT3 , 10X , 4KX C 1 ) , 3X , 

C $4KY ( 1 ) , 8X , 4HX (2 ) , BX , 4HY (2 ) , SX , 4HX (3 ) , 8X , 4L (3)) 

17 FORMAT (4X, 414, 8F10. 5) 

DO 10 J=l, LMENTS 

WRITE (10 , 11) J, (JT(3000*(I-l)-fJ) ,1=1,4) , (X(I,J) , Y(I , J) ,1=1 ,4) 
11 F0RMAT(1X,5I5,3X,SF12.4) 

10 CONTINUE 

READ (IN, 300) IBO 
300 FORMAT (50X, 14, 46X) 

READ(IN,301) (NBO(I) ,1=1, IBO) 

C WRITE (10, 302) 

DO 333 1=1, LMENTS 

READ (IN, 444) (NOO(I , J) , J=1 , 7 ) 

333 CONTINUE 

C WRITE (10, £03) (NBO (I ), 1=1, IBO) 

302 FORMATC//, ' *** BOUNDARY NUMBERS ***' ) 

CALL SETUP 
NTBAN=IDIFF+1 
WRITE (10, 202)NTBAN 

202 FORMAT (// , 2X,25HTHE ORIGINAL BANDWIDTH IS, 14,//) 

CALL OPTNUM 

IF ( ID IFF . LE . MINMAX ) GO TO 115 
M I B AN=MINM AX+1 


A 39 


ORBGTOL PAGE G£f 
OF POOR QUALITY 


WRITE(IQ, 198) 

198 FORMAT ( 1H1/ / / , IX , 5 C3HOLD ,3k, 3HNEW , 7X ) , / , IX , 

§ 5 (4HN0DE , 2X , 4HN0DE , 6X) ) 

WRITE (7, 200) (J, JNT(J) ,J=1, NODES) 

200 FORMAT (5(15, IX ,15,5X3) 

C WRITE (10, 13) 

DO 180 J=1 , LMENTS 
DO 555 11=1,4 
I I 1=11+2 

555 N00( J, III)=JNT(JT(300Q*(II-1)+J) ) 

C 

180 CONTINUE 

C WRITE(IO, 11) J, (JNT(JT(3000*(I-1)+J) ) , 1=1,4) , 

C $(X(I,J) , Y (I, J), 1=1, 4) 

WRITE ( I 0 , 20 1 ) MIB AN 

201 FORMAT ( / / , 2X , * THE NEW BANDWIDTH IS ' ,13,/) 

C WRITE (IP , 310 ) NODES , LMENTS, MIBAN 

C DO 181 J=l, LMENTS 

C81 WRITE (IP, 18) J, (JNT(JT(30Q0*(I-1)+J)) , 1=1,4) , (X(I, J) , 
C $Y(I, J) ,1=1,4) 

C8 F0RMAT(5I5 , 8F10 .4) 

C WRITE (10, 302) 

C WRITE(XP, 300) IBO 

C WRITE (IP ,301) (JNT(NBO(I) ) , 1=1 , IBO) 

C WRITE(IO, 303) ( JNT(NBO(I) ) , 1=1 , IBO) 

303 FORMAT (15 14) 

GO TO 117 
115 CONTINUE 

WRITE (10, 101) 

117 CONTINUE 

101 F0RMAT(///,2X, 'THE ORIGINAL BANDWIDTH IS MINIMUM 1 ) 

WRITE ( I P , 3 1 0 ) NODES , LMENTS , NTB AN 
310 FORMAT (3 15) 

C DO 190 J=l, LMENTS 

C190 WRITE (IP, 18) J, (JT(3000*(I-1)+J) , 1=1,4) , 

C 1(X(I,J) ,Y(I,J) ,1=1,4) 

C WRITE (IP, 300) IBO 

C WRIT£(IF, 301) (NBO(I), 1=1, IBO) 

301 FORMAT (2015) 

C- C! 

444 FORMAT (6 15, A4) 

LLL = N00(1 , 7) 

C 

L0L=0 

MMM=0 

NNN=1 

DO 666 11=1, LMENTS 

IF ( LLL.EQ.NOO(H j 7) )GO TO 1888 

NNN=NNN-1 

L0L=L0L+NNN 

MMM=MMM+1 

NTT (MMM)=LOL 

C WRITE (10, 888) NNN,NTT(MMM) 

LLL=N00(II , 7) 


o n o 
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NNN=1 

1888 CONTINUE 

WRITE (10, 777 )NNN, (NOO(II, J) , J=l,7) , (X(I,II) ,Y(I, II) , 1=1,4) 
C WRITE (7, 1012)NQ0(II,2) ,N00(II,6) 

WRITE(7 , 199) (NOO(IIjJ) , J=2,6) , (X(I,II),Y(I,II) ,1=1,4) 

IF (II. NE . LMENTS) GO TO 1889 

LOIr=LOL+NNN 

MMM=MMM+1 

NTT(MMM)=LOL 

C WRITE (10, 888) NNN,NTT(MMM) 

1889 NNN=NNN+1 
666 CONTINUE 

C 

C1012 FORMAT (15 * A5 , 15 , A5 , 15 , A5 , 15 , A5 , 15 , A5 , 15 , A5 , 15 , A5 ) 

777 FORMAT (2X, 715 , IX, A8 , 8F10 .5) 

199 FORMAT (5I4,8F7.4) 

888 FORMAT (// * NUMBER OF ELEMENTS IN THIS REGION =*,2110/) 

400 CONTINUE 


STOP 

END 
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TRASFER PROGRAM 


55 

44 


C 


c 

5 


1 


1 

1 

1 

1 

1 

1 

1 

1 


DIMENSION XC4S2) ,Y(482) ,T(482) ,BETA(376) , PROB (3, 7) , ALFA (3 ,3) , 
£1(150) ,TA(482) ,TG(482) 

INTEGER NM,I,J,K, II ,IE(376,5) , III , IT,MID(376) ,IG(482) , IS(482) 
CHARACTERS TIP£(376) ,PT(150) 

DOUBLE PRECISION PROB , ALFA , BETA 
DATA PROB /770000. 0,97500. 0,1000.0, 

164000.0,132600.0,1500.0, 

1330000.0,1950000.0,3000.0, 

0.440,0.287,0.40, 

0.166,0.022,0.02, 

0.255,0,02,0.02, 

76000.0,14000.0,500.0/ 

DATA ALFA /0. 00002705 , 0 . 00003490, 0 . 0000005 , 

0 . 00003036 , 0 . 0000259,0. 0000005 , 
0.00001457,0.0000021,0.0000005/ 

DATA BIRT/'BIRT'/.BIRL/'BIHL'/, FI/' FI'/, PAP /'PAP'/ 


DO 3 I = 1,3 

WRITE(7 ,55) (PROB (I, J) , J « 1,7) 
WRITE (7, 44) (ALFA(I , J) , J =1,3) 
CONTINUE 

FORMAT (3F 10 . 1 , 3F10 . 8 , F10 . 1) 
FORMAT (3F 10. 8) 

PI = 3.1415927 


TI = 70.0 


DO 5 I as 1,389 
IF(I.LE.140)TIPE(I) ='BIR' 
IF((I .GE. 141) .AND. (I .LE . 151) ) 
IF( (I .GE. 152) .AND. (I.LE.175)) 
IFfCI.GE. 176). AND. (I.LE.319)) 
IF((I . GE . 320) .AND . (I .LE. 343) ) 
IF((I . GE. 344) .AND. (I.LE.376)) 
IF((I.GE.377) .AND. (I.LE.389)) 
CONTINUE 


TIPE(I) ='FIR' 
TIPE(I) = 1 PAP * 
TIPE(I) ='FIR 1 
TIPE(I) =’PAP ! 
TIPE(I) = , FIR’ 
TIPE(I) =' BIR 1 


DO 51 KK « 1,389 

IF(KK.GE.l. AND.KK.LE.il) BETA (KK)=- . 2618 
IF(KK.GE. 12. AND .KK.LE . 17) BETA(KK)=- . 15707 
IF (KK . GE . 18 . AMD . KK . LE . 29 ) BETA(KK)=0.0 
IF(KK.GE.30.AND.KK.LE.33) BETA(KK)=. 38394 
IF (KK . GE . 34 . AND . KK . LE . 3 7 ) BETA (KK)=. 17452 
IF (EK.GE . 38. AND .KK.LE .41) BETA (KK)=1. 22164 
IF (KK . GE . 42 . AND . KK . LE . 45 ) BETA(KK)=.69808 
IF (KK . GE . 46 . AND . KK . LE . 49 ) BETA (KK}=2. 4433 
IF (KK . GE . 5 0 . AND . KK . LE . 5 3 ) BETA (KK)=1 . 91972 
IF{KK-GE. 54. AND. KK.LE. 57) BETA(KK)=2 . 9668 
IFCKK.GE. 58. AND, KK.LE. 61) BETA (KK) =2. 7574 
IF(KK,GE. 62.AND. KK.LE. 73) BETA (KK)=3 . 14136 
IFCKK.GE. 74. AND. KK.LE. 79) BETA (KK)-3 .29843 
IFCKK.GE. 80. AND. KK.LE. 101) BETA (KK)=3 . 40314 
IFCKK.GE. 102. AND. KK.LE. 107) BETA (,KK)=3 . 29843 
IFCKK.GE. 108. AND. KK.LE. 115) BETA (KK) =3 . 14136 
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51 


C 

200 


10 

777 


12 

666 


BETA(KK)=0.0 
BETA(KK)=-. 15707 
BETA(KK)=“.2618 
BETA(KK)=~.15707 
BETA (KK) =0.0 
BETA(KK)=. 38394 
BETA(KR)=. 17452 
BETA(KK)=. 38394 
BETA(KK)=. 17452 
BETA (KK)=1. 2216 
BETA (KK)=. 69808 
BETA (KK) =1.2216 
BETA (KK)=. 69P-08 
BETA (KK)=2. 4433 
BETA (KK)=1. 91972 
BETA(KK)=2 . 4433 
BETA(KK)=1 .91972 
BETA (KK)=2. 9668 
BETA (KK)=2. 7574 
BETA (KK)=2. 9666 
BETA(KK)=2 . 7574 
BETA (KK)=3. 14136 
BETA (KK)=3. 29843 
BETA (KK)=3 . 403 14 
BETA (KK) =0.0 
BETA(KK)=1. 57079 


25 

50 


IFCKK.GE.116.AND.KK.LE.123) 

IFCKK.GE. 124.AND.KK.LE. 129) 

IF(KK.GE. 130.AND.KK.LE, 163) 

IFCKK.GE.164.AND.KK.LE. 175) 

IFCKK.GE.176.AND.KK.LE. 215) 

XFCKK.GE.216.AND.KK.LE. 219) 

IF (KK . GE . 220 . AND . KK . BE . 223 ) 

IF (KK . G£ . 224 . AND . KK . LE . 2 2 7 ) 

IF CKK.GE.228.AND.KK.LE. 231) 

IF (KK . GE . 232 . AND . KK. LE . 235 ) 

IFCKK.GE. 236. AND.KK.BE. 239) 

IF (KK . GE . 240 . AND .KK . LE . 243 ) 

IFCKK . GE . 244 .AND .KK . LE . 247 ) 

IFCKK.GE. 248. AND.KK.BE. 251) 

IF CKK.GE.252.AND.KK.LE. 255) 

IF(KK. GE.256.AND.KK.LE. 259) 

IF CKK.GE.260.AND.KK.LE. 263) 

IF(KK . GE . 264 . AND . KK. BE . 267 ) 

IF (KK , GE . 268 . AND . KK . LE . 27 1 ) 

IF CKK.GE.272.AND.KK.BE. 275) 

IF CKK.GE.276.AND.KK.LE. 279) 

IF CKK.GE.2S0.AND.KK.BE. 319) 

IF CKK.GE.320.AND.KK.LE. 331) 

IF CKK.GE.332.AND.KK.LE. 354) 

IF CKK.GE.355.AND.KK.LE. 376) 

IF CKK.GE.377.AND.KK.LE. 389) 

CONTINUE 

READ(5 ,20Q) (IG(I) , 1=1,482) 

WRITE (6, 200) (IG(I) , 1=1,482) 

READ (5 ,200) (IS(I) , 1=1,482) 

W r RITE(6 , 200) (IS(I) ,1=1,482) 

FORMAT (5 (6X , 15 , 5X) ) 

DO 10 L = 1,376 

READ(5, 777) NM,I, J,K,N,X(I) ,Y(I) ,X(J) , Y(J) ,X(K) , Y(K) ,X(N) ,Y(N) 
WRITE (6, 777) NM, I, J,K,N,X(I) , Y(I) ,X(J) ,Y(J) ,X(K) ,Y(K) ,X(N) ,Y(N) 
IECL,1) = NM 
IE(L,2) = I 

IE(L,3) = J 
IE(L,4) = K 
IE(L,5) = N 
CONTINUE 

FORMAT (5 14, 8F7. 4) 

READ (5 , 666) (T(X), I = 1,482) 

WRITE (6,666) (T(I), I = 1,482) 

DO 12 I = 1,482 
TG(IGCI) )=T(IS (I) ) 

FORMAT (9F8. 4) 

MID Cl) =1 
DO 50 I = 2,376 

IF((TIPE(I) -NE .TIPE (1-1) ) .OR. (BETA (I ) .NE. BETA(I-l) ) )GQ TO 25 
MID (I) = MID ( I - 1 ) 

GO TO 50 

MID (I) = MID (1-1) -f- 1 
CONTINUE 
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IT = 0 

DO 70 I = 1,376 
IF (MID (I) .GE. IT) GO TO 15 

15 IT = MID (I) 

BI(IT) = BETA(I) 

PT(IT) = TIPE(I) 

WRITE (6 , 16)1 , TIPE(I) , BETA (I), IT 
C WRITE (6, 66) IT,PT(IT) 

70 CONTINUE 

16 FORMAT(I5 ,5X,A4,F10 .4, 15) 

DO 30 I * 1,482 

TA(I) = 90.0 + TG(I) 

WRITE (7, 555) I ,X(I) ,Y(I) ,TA(I) 

30 CONTINUE 

555 FORMAT (15 ,40X,2F10.4,5X,F10.4) 

DO 80 1=1, IT 

C WRITE (6, 66) I, PT(I) 

BII = BI (I)*180 .0/PI 
WRITE (7, 11) I, BII 
IF(PT(I).EQ. 'BIR') GO TO 7 
IF(PT(I) .EQ. *FIR f ) GO TO 17 
IF(PT(I).EQ. 'PAP') GO TO 27 
7 WRITE (7, 22) TI , (PR0B(1, J) , J=l,7) 

WRITE (7 ,33) (ALFA(1, J) , J =1,3) 

GO TO 80 

17 WRITE(7 ,22) TI , (PROB(2, J) , J=l,7) 

WRITE (7 ,33) (ALFA(2, J) , J =1,3) 

GO TO SO 

27 WRITE (7, 22) TI , (PROB(3, J) , J=l,7) 

WRITE (7 , 33) (ALFA(3, J) , J =1,3) 

80 CONTINUE 

66 FORMAT (15 , 5X, A4) 

11 F0RMAT(I5,25X,F10.4) 

22 FORMAT(4F10 . 1 ,3F10. 8,F10 . 1) 

33 FORMAT (3F 10 .8) 

TH=1. 000 

N= 4 

TR = 70.00 
C M =1 

DO 40 I =1,376 
WRITE (7, 444) (IE(I , J) , J=1 ,5) ,MID(I) ,TR, N, TH 
40 CONTINUE 

444 FORMAT (6 15 ,F10 .4, 10X, 15 ,5X,F10.4) 

STOP 

END 
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//GQ.SYSIN DD * 

THERMAL STRESSES IN THE BLADE 

USING PLANE STRAIN APPROACH 


496 1 1 

1 -1 

-1 -1 

“1 4.3519 

13.0728 

147.8842 

2 


3.8200 

12.5000 

144.3767 

3 


3.2448 

14.0923 

150.2360 

4 


2.6075 

13.3675 

147.7226 

5 


5.0275 

13.6888 

150.4952 

6 


4.0394 

14.8344 

152.0367 

7 


3.4319 

11.9703 

139.9296 

8 


2.1273 

12.6598 

144.3721 

9 


2.1441 

15.1056 

154.4340 

10 


1.4000 

14.2300 

154.0150 

11 


3.0588 

15.9725 

154.7096 

12 


0.8266 

13.3456 

153.5614 

13 


5.8469 

14.3478 

152.3099 

14 


4.9911 

15.5936 

153.2852 

15 


4.1441 

16.8306 

154.9098 

16 


3. 1875 

11.4838 

134.7187 

17 


1.8044 

11.9694 

139.7209 

18 


0.4237 

12.4525 

152.5827 

19 


2.0344 

15.2103 

154.7712 

20 


1.2800 

14.3200 

154.5201 

21 


2.9625 

16.0887 

154.9232 

22 


0.6994 

13.4178 

154.2866 

23 


4.0644 

16.9553 

155.0355 

24 


0.2925 

12,5037 

153.5707 

25 


6.8100 

15.0500 

153.5261 

26 


6.1000 

16.3700 

154.1107 

27 


5.4000 

17.6800 

155.0395 

28 


5.3400 

17.8100 

155.1139 

29 


3.0869 

11.0403 

129.1960 

30 


1.6386 

11.2961 

133.0223 

31 


0. 1916 

11.5506 

150.4510 

32 


0.0594 

11.5778 

154.5046 

33 


7.2091 

15.3441 

153.9100 

34 


6.5538 

16.6739 

154.3554 

35 


5.9059 

18.0200 

155.0804 

36 


5.8544 

18.1544 

155.1307 

37 


3.1300 

10.6400 

124.1252 

38 


1.6300 

10.6400 

123.5315 

39 


0.1300 

10.6400 

123.5111 

40 


0.0000 

10.6400 

123.9204 

41 


7.8038 

15.6287 

154.2687 

42 


7.2100 

16.9706 

154.5677 

43 


6.6213 

18.3550 

155.1068 

44 


6.5775 

18.4925 

155.1460 

45 


3,0397 

10.2559 

119,1247 

46 


1.5559 

10.0748 

115.2346 

47 


0.0772 

9.8925 

97.1668 


A 4^ 


48 


- 0,0594 

9.8781 

92.0556 

49 


8.5941 

15.9040 

154.5623 

50 


8.0688 

17.2601 

154.7453 

51 


7.5459 

18.6850 

155.1289 

52 

ORIGINAL mm 

7.5094 

18.8243 

155.1588 

53 

OF POOR QUALITY 

3.0913 

9.8112 

113.2342 

54 


1.6412 

9.4569 

108. 1028 

55 


0.2013 

9.1000 

93.2964 

56 


0.0625 

9.0725 

91.9792 

57 


9.5800 

16.1700 

154.7605 

58 


9.1300 

17.5425 

154.8743 

59 


8.6800 

19.0100 

155.1474 

60 


8.6500 

19.1500 

155.1694 

61 


3.2847 

9.3059 

107.4809 

6 ? 


1.8859 

8.7861 

102.8803 

63 


0.5022 

8.2625 

91.9885 

64 


0.3656 

8.2231 

91.1432 

65 


10.7616 

16.4266 

154.8534 

66 


10.3938 

17.8177 

154.9407 

67 


10.0234 

19.3300 

155.1605 

68 


9.9994 

19.4693 

155.1775 

69 


3.6200 

8.7400 

102.4331 

70 


2.2900 

8.0625 

98.9507 

71 


0.9800 

7.3800 

91.4340 

72 


0.8500 

7.3300 

90.8650 

73 


12.1388 

16.6737 

154.8177 

74 


11.8600 

18.0856 

154.9131 

75 


11.5763 

19.6450 

155.1579 

76 


11.5575 

19.7825 

155.1759 

77 


4.0972 

8.1134 

98.3440 

78 


2.8534 

7.2861 

95.8928 

79 


1.6347 

6.4525 

90.9115 

80 


1.5156 

6.3931 

90.5528 

81 


13.7116 

16.9116 

154.5817 

82 


13.5288 

18.3464 

154.7357 

83 


13.3384 

19.9550 

155 . 1338 , 

84 


13.3244 

20.0894 

155.1635 

85 


4.7163 

7.4263 

95.3039 

86 


3.5763 

6.4569 

93.6887 

87 


2.4663 

5.4800 

90.5535 

88 


2.3625 

5.4125 

90.3444 

89 


15.4800 

17.1400 

153.8904 

90 


15.4000 

18.6000 

154.2581 

91 


15.3100 

20.2600 

155.0557 

92 


15.3000 

20.3900 

155.1067 

93 


5,4772 

6.6784 

93-2149 

94 


4.4584 

5.5748 

92.1772 

95 


3.4747 

4.4625 

90.3037 

96 


3.3906 

4.3881 

90.1936 

97 


16.3600 

17.2000 

153.0317 

98 


16.6300 

17.7400 

152.9341 

99 


15.8000 

17.8000 

153.7650 

100 


15.8025 

18.9350 

154.2794 

101 


15.8000 

20.2600 

155 . 0328 


A4-C? 
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102 

15.8000 

20.3900 

155.0989 

103 

6.3800 

5.8700 

91.8505 

104 

5.5000 

4.6400 

91.1968 

105 

4.6600 

3.4000 

90.1492 

106 

4.6000 

3.3200 

90.0924 

107 

16.8200 

17.2300 

152.0987 

108 

16.9700 

18.0400 

152.8914 

109 

17.0800 

18.8300 

153.7508 

110 

16.8200 

18.2800 

153.2712 

111 

16.3800 

18.3600 

153.6377 

112 

16.3400 

19.2200 

154.3084 

113 

16.3000 

20.2600 

155.0112 

114. 

16.3000 

20.3900 

155.0851 

115 

6.7350 

5.5134 

91.4067 

116 

5.9581 

4.2261 

90.8833 

117 

5.2137 

2.9275 

90.1065 

118 

5.1578 

2.8425 

90.0704 

119 

17.2200 

17.2200 

150.9112 

120 

17.2200 

18.3000 

153.1369 

121 

17.2200 

19.3400 

154.2458 

122 

16.7600 

19.0000 

154.0092 

123 

17.2200 

19.4600 

154.3388 

124 

16.7600 

19.5750 

154.4896 

125 

16.7600 

20.2600 

154.9978 

126 

16.7600 

20.3900 

155.0772 

127 

7.2900 

5 . 1663 

90.9979 

128 

6.6150 

3.8294 

90.6383 

129 

5.9650 

2.4800 

90.0829 

130 

5.9137 

2.3900 

90.0547 

131 

17.8400 

17.5400 

151.4541 

132 

17.9050 

18.4900 

153.2845 

133 

17.2200 

15.6538 

144.1907 

134 

17.8725 

15.7737 

144.6731 

135 

17.9700 

19.3600 

154.2331 

136 

17.9700 

19.5000 

154.3391 

137 

17 . 2200 

19.8400 

154.6544 

138 

17.9700 

19.8750 

154.6646 

139 

17.2200 

20.2600 

154.9879 

140 

17.2200 

20.3900 

155 . 0715 

141 

8.0450 

4.8284 

90.6712 

142 

7.4706 

3.4498 

90.4450 

143 

6.9137 

2.0575 

90.0612 

144 

6.8678 

1.9625 

90.0401 

145 

18.4700 

17.8600 

152.5973 

146 

18.4700 

18.6700 

153.6167 

147 

18.4700 

15.8938 

145.1270 

148 

18.4700 

19.3600 

154.2 S 30 

149 

17.2200 

14.1600 

138.0302 

150 

17.9050 

14.1600 

138.0316 

151 

18.4700 

14.1600 

138.0318 

152 

18.4700 

19.5400 

154.4113 

153 

18.4700 

19.9100 

154.7148 

154 

17.9700 

20.2600 

154.9807 

155 

18.4700 

20.2600 

154.9871 


156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 
209 



17.9700 

20.3900 

155.0665 


9.0000 

4.5000 

90.4415 


8.5250 

3.0875 

90.3062 

cmsGML pass m 

8.0600 

1.6600 

90.0439 

OF POOR QUALITY 

8.0200 

19.4700 

1.5600 

17.8700 

90.0285 

153.7898 


19.4700 

18.6550 

154.0946 


19.4700 

19.3700 

154.5283 


19.4700 

19.5200 

154.6021 


17.2200 

12.7388 

132.1555 


17.9375 

12.6987 

132.0249 


18.4700 

12.6587 

131.8744 


19.4700 

19.8925 

154.8245 


19.4700 

20.2600 

155.0414 


18.4700 

20.3900 

155.0680 


19.4700 

20.3900 

155.1040 


10.1550 

4.1809 

90.3129 


9.7781 

2 . 7423 

90.2227 


9.4038 

1.2875 

90.0314 


9.3703 

1.1825 

90.0198 


20.4700 

17.8800 

154.1075 


20.4700 

18.6400 

154.3112 


20.4700 

19.3800 

154.7841 


20.4700 

19.5100 

154.8236 


20.4700 

19.8800 

154.9594 


17.2200 

11.3900 

126,2683 


17.9700 

11.3900 

126.7595 


18.4700 

11.3900 

126.8414 


20.4700 

20.2600 

155.0976 


20.4700 

20.3900 

155.1358 


11.5100 

3.8713 

90.2962 


11.23Q0 

2.4144 

90.2166 


10.9450 

0.9400 

90.0313 


10.9187 

0.8300 

90.0196 


21.7062 

19.3500 

155.1174 


21.7012 

19.4887 

155.1178, 


21.7181 

19.8631 

155.1418 


21.7325 

20.2350 

155.1785 


17.2200 

10.6400 

124.3878 


17.9700 

10.6400 

124.4096 


16.4700 

10.6400 

124.3720 


16.4700 

11.3900 

125.3677 


18.4700 

10.6400 

124.4140 


21.7387 

20.3550 

155.1893 


13.0650 

3.5709 

90.4361 


12.8806 

2.1036 

90.3264 


12.6838 

0.6175 

90.0457 


12.6653 

0.5025 

90.0261 


22.9400 

19.3000 

155.1814 


22.9400 

19.4400 

155.1815 


22.9675 

19.8175 

155.1866 


22.9900 

20. 1800 

155.1946 


23.0000 

20.3000 

155.1970 


17.2200 

9.8500 

122.3997 
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210 

17.9700 

9,8500 

121.9564 

211 

16.4700 

9.8500 

123.2855 

212 

18.4700 

9.8500 

121.8808 

213 

14.8200 

3.2800 

90.9077 

214 

14.7300 

1.8100 

90.6688 

215 

14.6200 

0.3200 

90.1157 

216 

14.6100 

0.2000 

90.0768 

217 

24.1712 

19.2300 

155.1958 

218 

24.1862 

19.3637 

155.1959 

219 

24.2181 

19.7431 

155.1970 

220 

24.2425 

20.0950 

155.1988 

221 

24.2537 

20,2250 

155.1994 

222 

17.2200 

8.1750 

115.1947 

223 

17.9400 

8.2131 

115.3302 

224 

18.4700 

8.2537 

115.4882 

225 

15.2600 

2.7000 

91.1052 

226 

15.2650 

1.4925 

90.7231 

227 

16.0200 

3.2000 

91.7380 

228 

16.3000 

2.6500 

91.8801 

229 

15.2600 

0.2800 

90.1457 

230 

15.2400 

0.1500 

90.0878 

231 

25 . 4000 

19 . 1400 

155.1987 

232 

25.4400 

19.2600 

155.1988 

233 

25.4700 

19.6400 

155.1991 

234 

25.4900 

19.9800 

155.1997 

235 

25.5000 

20.1300 

155.1999 

236 

17.2200 

6.5000 

108.2843 

237 

17.9100 

6.4975 

108.2734 

238 

18.4700 

6.5000 

108.2838 

239 

15.9800 

2.1000 

91.3425 

240 

15.9400 

1.1600 

90.7478 

241 

16.6000 

2.1200 

91.7607 

242 

15.9000 

0.2300 

90.1710 

243 

16.6000 

3.1600 

92.6067 

244 

16.7350 

2.3650 

92.1014 

245 

16.8800 

1.5800 

91.3814 

246 

15.8800 

0.1000 

90.1029 

247 

27.4722 

18.9144 

155.1998 

248 

27.5022 

19.0455 

155.1998 

249 

27.5397 

19.4275 

155.1999 

250 

27.5722 

19.7844 

155.2000 

251 

27.5855 

19.8989 

155,2000 

252 

17.2200 

4.8250 

101.3809 

253 

17.8800 

4.7031 

100.8879 

254 

18.4700 

4.5887 

100.4528 

255 

16.6000 

1.5100 

91.2262 

256 

16.5800 

0.8375 

90.6678 

257 

16.5600 

0.2000 

90.1982 

258 

16.5500 

0.0700 

90.1203 

259 

17.2200 

3.1500 

94.2286 

260 

17.2200 

2.1000 

92.0724 

261 

17.2200 

1.0500 

90.9383 

262 

17.2200 

0.9300 

90.8391 

263 

29.5188 

18.6544 

155.2000 


ft ^ 
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264 


29.5422 

18,7922 

155.2000 

265 


29.5856 

19.1733 

155,2000 

266 


29.6289 

19.5444 

155.2000 

267 

268 
269 

original pags m 
OF POOR QUALITY 

29.6455 

17.8500 

18.4700 

19.6355 

2.8300 

2.5200 

155.2000 

93.6769 

92.5586 

270 


17.2200 

0.5200 

90.4926 

271 


17.2200 

0,1700 

90.2098 

Til 


17.2200 

0.0500 

90.1317 

273 


17.9100 

1.9350 

91.9496 

274 


17.9700 

1.0300 

90,9540 

275 


17.9700 

0,9100 

90.8619 

276 


17.9700 

0.5200 

90.5194 

277 


31.5400 

18.3600 

155.2000 

278 


31,5600 

18.5000 

155.2000 

279 


31.6075 

18.8775 

155.2000 

280 


31.6600 

19.2600 

155.2000 

281 


31.6800 

19 . 3400 

155,2000 

282 


18,4700 

1 . 7800 

91.6346 

283 


19.4700 

2.5400 

91.3945 

284 


19.4700 

1.7900 

91.1279 

285 


17.9700 

0.1600 

90.2209 

286 


17.9700 

0.0300 

90.1343 

287 


18.4700 

1 . 0200 

90.8996 

288 


18.4700 

0.9000 

90.8143 

289 


18.4700 

0.5200 

90.5006 

290 


18.4700 

0.1400 

90.2044 

291 


33,5355 

18.0311 

155.2000 

292 


33.5555 

18.1689 

155.2000 

293 


33.6055 

18.5400 

155.2000 

294 


33.6655 

18.9311 

155.2000 

295 


33.6888 

19.0122 

155.2000 

296 


19.4700 

1.0300 

90.6696 

297 


20.4700 

2.5600 

91.1005 

298 


20.4700 

1.8000 

90.9125 

299 


20.4700 

1.0400 

90.4442. 

300 


18.4700 

0.0200 

90.1300 

301 

' 

19.4700 

0.9000 

90.6054 

302 


19.4350 

0.5150 

90.3788 

303 


19.4100 

0.1400 

90.1577 

304 


19.4100 

0.0200 

90.0994 

305 


35.5055 

17.6678 

155.2000 

306 


35.5289 

17.7989 

155.2000 

307 


35.5797 

18.1608 

155.2000 

308 


35 . 6455 

18.5578 

155.2000 

309 


35.6722 

18.6522 

155.2000 

310 


20.4700 

0.9100 

90.3968 

311 


21.2187 

1.0387 

90.1452 

312 


21.2175 

0.9162 

90.1430 

313 


20.4700 

0.5200 

90.2397 

314 


20.4900 

0.1500 

90.0937 

315 


20.4700 

0.0400 

90.0601 

316 


37.4500 

17.2700 

155.2000 

317 


37.4800 

17.3900 

155.2000 
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318 


37.5299 

17.7400 

155,2000 

319 


37.6000 

18.1400 

155.2000 

320 


37.6300 

18.2600 

155.2000 

321 

pags m 

21.2100 

0.5300 

90.1030 

322 

OF POOR QUALITY 

22.2900 

1.0800 

90.0357 

323 

22.2900 

0.9600 

90.0355 

324 


22.2850 

0.5800 

90.0275 

325 


21.2175 

0.1538 

90.0439 

326 


21.2100 

0.0375 

90.0284 

327 


37.8422 

17.2166 

155.2000 

328 


37.8555 

17.3455 

155.2000 

329 


37.8680 

17.7036 

155.2000 

330 


37.8988 

18.1066 

155.2000 

331 


37.9122 

18.2377 

155,2000 

332 


22,2900 

0.2000 

90.0110 

333 


23.6837 

1,1637 

90.0063 

334 


23.6875 

1.0413 

90.0062 

335 


23.6950 

0.6700 

90,0048 

336 


23.7075 

0 . 2887 

90,0019 

337 


22.2900 

0.0800 

90.0062 

338 


38,8589 

16.9933 

155.2000 

339 


38.8622 

17.1289 

155,2000 

340 


38.8555 

17.4911 

155.2000 

341 


38.8655 

17.8933 

155.2000 

342 


38.8689 

18.0311 

155.2000 

343 


23.7100 

0.1675 

90.0010 ’ 

344 


25.4000 

1.2900 

90.0013 

345 


25.4100 

1.1600 

90.0013 

346 


25.4400 

0.8000 

90.0009 

347 


25.4700 

0.4200 

90.0002 

348 


25 .4700 

0.3000 

90.0001 

349 


40.5000 

16.6000 

155.2000 

350 


40.5000 

16.7400 

155.2000 

351 


40.4925 

17.1025 

155.2000 

352 


40.5000 

17.5000 

155.2000 

353 


40.5000 

17.6400 

155.2000 

354 


27,3933 

1.5367 

90.0002 

355 


27.4111 

1.4111 

90.0002 

356 


27.4577 

1.0478 

90.0001 

357 


27.5011 

0.6711 

90.0000 

358 


27.5155 

0.5456 

90.0000 

359 


42.7655 

16.0367 

155.2000 

360 


42.7689 

16.1789 

155.2000 

361 


42.7788 

16.5378 

155.2000 

362 


42.8022 

16.9267 

155.2000 

363 


42.8055 

17.0644 

155.2000 

364 


29.3733 

1.8133 

90.0000 

365 


29.3978 

1.6911 

90.0000 

366 


29.4577 

1.3244 

90.0000 

367 


29.5111 

0.9511 

90.0000 

368 


29.5355 

0.8222 

90.0000 

369 


45 . 6555 

15.3033 

155.2000 

370 


45.6689 

15.4456 

155.2000 

371 

C ~ 2^ 

45.7147 

15.7970 

155.2000 
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372 

45,7722 

16.1733 

155.2000 

373 

45.7855 

16.3044 

155.2000 

374 

31.3400 

2.1200 

90.0000 

375 

31.3700 

2.0000 

90.0000 

376 

31.4400 

1.6300 

90.0000 

377 

31.5000 

1.2600 

90.0000 

378 

31.5300 

1.1300 

90.0000 

379 

49.1700 

14.4000 

155.2000 

380 

49 . 2000 

14.5400 

155.2000 

381 

49.2999 

14.8800 

155.2000 

382 

49 .4100 

15.2400 

155.2000 

383 

49.4400 

15.3600 

155.2000 

384 

33.2933 

2.4567 

90.0000 

385 

33.3277 

2.3378 

90.0000 

386 

33.4044 

1.9644 

90.0000 

387 

33.4677 

1.5978 

90.0000 

388 

33.4989 

1.4689 

90,0000 

389 

50.4525 

14.0200 

155.2000 

390 

50.5025 

14.1600 

155.2000 

391 

50.6056 

14.5037 

155.2000 

392 

50.7162 

14.8625 

155.2000 

393 

50.7450 

14.9675 

155.2000 

394 

35.2333 

2.8233 

90.0000 

395 

35.2711 

2.7044 

90.0000 

396 

35.3511 

2.3278 

90.0000 

397 

35.4144 

1,9644 

90.0000 

398 

35.4422 

1.8389 
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238 

224 

223 

44 

70.0000 

4 

1.0000 

359 

252 

253 

237 

236 

44 

70.0000 

4 

1.0000 

360 

253 

254 

238 

237 

44 

70.0000 

4 

1.0000 

361 

259 

268 

253 

252 

44 

70.0000 

4 

1.0000 

362 

268 

269 

254 

253 

44 

70.0000 

4 

1.0000 

363 

194 

195 

182 

181 

44 

70.0000 

4 

1.0000 

364 

195 

198 

183 

182 

44 

70.0000 

4 

1.0000 

365 

209 

210 

195 

194 

44 

70.0000 

4 

1.0000 

366 

210 

212 

198 

195 

44 

70.0000 

4 

1.0000 

367 

196 

194 

181 

197 

44 

70.0000 

4 

1.0000 

368 

211 

209 

194 

196 

44 

70.0000 

4 

1.0000 

369 

133 

134 

131 

119 

44 

70.0000 

4 

1.0000 

370 

134 

147 

145 

131 

44 

70.0000 

4 

1.0000 

371 

149 

150 

134 

133 

44 

70.0000 

4 

1.0000 

372 

150 

151 

147 

134 

44 

70.0000 

4 

1.0000 

373 

165 

166 

150 

149 

44 

70,0000 

4 

1.0000 

374 

166 

167 

151 

150 

44 

70.0000 

4 

1,0000 

375 

181 

182 

166 

165 

44 

70.0000 

4 

1.0000 

376 

182 

183 

167 

166 

44 

70.0000 

4 

1 .0000 

377 

434 

485 

484 

435 

45 

70.0000 

4 

1.0000 

378 

433 

486 

485 

434 

45 

70.0000 

4 

1.0000 

379 

430 

487 

486 

433 

45 

70.0000 

4 

1.0000 

380 

429 

488 

487 

430 

45 

70.0000 

4 

1.0000 

381 

432 

489 

488 

429 

45 

70.0000 

4 

1.0000 

382 

483 

490 

489 

432 

45 

70.0000 

4 

1.0000 

383 

478 

491 

490 

483 

45 

70.0000 

4 

1.0000 


j. 

l! 


384 

476 

492 

491 

478 

45 

70.0000 

4 

1.0000 

385 

477 

493 

492 

476 

45 

70.0000 

4 

1.0000 

386 

480 

494 

493 

477 

45 

70.0000 

4 

1.0000 

387 

481 

495 

494 

480 

45 

70.0000 

4 

1.0000 

388 

482 

496 

495 

481 

45 

70.0000 

4 

1.0000 


1.0 
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RADIAL-TANGENTIAL STRESS PRGRAM 

DIMENSION SZETA(90) ,SEATA(90) ,SZT(90) ,BETA(90) 

REAL Sll, S22, S12, AA, BB 
INTEGER NN, NEL 
DO 20 KK - 1,90 

IF(KK.GE.l. AND. KK.LE. 11) BETA(KK)=- .2618 
IFCKK.GE. 12. AND. KK.LE. 17) BETA(KK)=- . 15707 
IF(KK.GE. 1S.AND.KK.LE.29) BETA(KK)=0.0 
IF(KK . GE . 30 . AND .KK . LE . 33) BETA (KK)=. 38394 
IFCKK.GE. 34. AND. KK.LE. 37) BETA(KK)=. 17452 
IFCKK.GE. 38. AND. KK.LE. 41) BETA (KK)=1 .22164 
IFCKK.GE. 42. AND, KK.LE. 45) BETA(KK)=.69808 
IF (KK . GE . 46 . AND . KK . LE , 49 ) BETA (KK)=2 . 4433 
IF (KK . GE . 50 . AND . KK. LE . 53) BETA(KIC)=1 .91972 
IFCKK.GE. 54. AND. KK.LE. 57) BETA (KK)=2. 9668 
IFCKK.GE. 58. AND. KK.LE. 61) BETA (KK)=2. 7574 
IFCKK.GE. 62. AND. KK.LE. 73) BETA(KK)=3. 14136 
IFCKK.GE. 74. AND. KK.LE, 79) BETA(KK)=3, 29843 
IFCKK.GE. 80. AND .KK.LE. 90) BETA (KK) -3.403 14 
20 CONTINUE 

READ(5, 11)NN 
11 FORMAT (15) 

WRITE (6,15) 

15 FORMAT (// , 10X , 3HNEL , 9X, 6HS -ZETA, 9X , 6HS-EATA, 9X , 10HTUE-ZET-ET, //) 

DO 10 1=1, NN 

READ (5, 22) NEL, Sll, S22, S12 

AA=(Sll+S22)/2 

BB= (Sll-S22)/2 

SZETA (NEL ) =AA+B B*COS (2*BETA (NEL) ) +S 12*S IN ( 2*BETA (NEL) ) 

SE ATA ( NEL ) =AA - B B*CO S (2*BETA(NEL) ) - S 1 2*S IN ( 2*BETA (NEL) ) 
SZT(NEL) =-BB*S IN (2*BETA (NEL) )+S12*C0S (2*BETA (NEL) ) 

WRITE ( 6 , 33 )NEL , SZETA (NEL) , SEATA(NEL) , SZT (NEL) 

10 CONTINUE 

22 FORMAT ( 15 , 3E 15 . 5 ) 

33 FORMAT (5X, 15 ,5X,3F 15 .3,/) 

STOP 

END 
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TEMPERATURE IN DEGREE F 


NODE NO. 

TEMP. 

NODE NO. 

TEMP. 

NODE NO. 

TEMP. 

NODE NO. 

TEMP. 

1 

64.771 

2 

64.015 

3 

64.434 

4 

64.520 

5 

64.923 

6 

63.561 

7 

54.372 

8 

57.723 

9 

64.710 

10 

60.236 

11 

64.287 

12 

65.035 

13 

62.583 

14 

49.721 

15 

44.719 

16 

49.930 

17 

54.377 

18 

64.910 

19 

62.037 

20 

57.884 

21 

63,571 

22 

65.114 

23 

60.451 

24 

43.022 

25 

39.196 

26 

65.039 

27 

63,285 

28 

64.111 

29 

60.495 

30 

64,505 

31 

65,080 

32 

65.131 

33 

33.511 

34 

33.532 

35 

34.125 

36 

64.355 

37 

62.310 

38 

63.526 

39 

63.910 

40 

33.920 

41 

65.107 

42 

65 . 146 

43 

7.167 

44 

25 . 235 

45 

29 . 125 

46 

64.568 

47 

64.269 

48 

2.056 

49 

65 . 129 

50 

65.159 

51 

3.296 

52 

18.103 

53 

23.234 

54 

64.745 

55 

64.562 

56 

1.979 

57 

65 . 147 

58 

65.169 

59 

1.989 

60 

1.143 

61 

12.880 

62 

17.481 

63 

64.874 

64 

64.760 

65 

65.161 

66 

65.177 

67 

1,434 

68 

0.865 

69 

8.951 

70 

12.433 

71 

64.941 

72 

64.853 

73 

65 . 158 

74 

65.176 

75 

0,911 

76 

0.553 

77 

5.893 

78 

8.344 

79 

64.913 

80 

64.818 

81 

65.134 

82 

65 . 163 

83 

0.554 

84 

0.344 

85 

3.689 

86 

5.304 

87 

64.736 

88 

64.582 

89 

65.056 

90 

65.107 

91 

64.258 

92 

0.304 

93 

0.194 

94 

2.177 

95 

3.215 

96 

63.890 

97 

65.033 

98 

65.099 

99 

64.279 

100 

63.765 

101 

0.149 

102 

1.197 

103 

0.092 

104 

1.851 

105 

63.032 

106 

65.011 

107 

65.085 

108 

64.308 

109 

63.638 

110 

62.934 

111 

0.107 

112 

0.883 

113 

1.407 

114 

0.070 

115 

62.099 

116 

64.998 

117 

65.077 

118 

64.490 

119 

64.009 

120 

63.271 

121 

62.891 

122 

0.083 

123 

0,638 

124 

0.998 

125 

0.055 

126 

60.911 

127 

64.988 

128 

65.071 

129 

64.654 

130 

64.339 

131 

63.751 

132 

64.246 

133 

63.137 

134 

0.061 

135 

0.445 

136 

0.671 

137 

Q.G4Q 

138 

61.454 

139 

54.191 

140 

54.673 

141 

64.981 

142 

65.067 

143 

64.665 

144 

64.339 

145 

64.233 

146 

63.285 

147 

0.044 

148 

0.306 

149 

0.441 

150 

0.029 

151 

62.597 

152 

55.127 

153 

48.030 

154 

48 .032 

155 

48.032 

156 

64.987 

157 

65 . 068 

158 

64.715 

159 

64.411 

160 

64.283 

161 

63.617 

162 

0.031 

163 

0.223 

164 

0.313 

165 

0.020 

166 

63.790 

167 

64.095 

168 

42.156 

169 

42 . 025 

170 

41,874 

171 

65.041 

172 

65.104 

173 

64.824 

174 

64.602 

175 

64.528 

176 

0.031 

177 

0.217 

178 

0.296 

179 

0.020 

180 

64.107 

181 

64.311 

182 

64.784 

183 

36.268 

184 

36.760 

185 

36.841 

186 

65.098 

187 

65.136 

188 

64.959 
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189 

64.824 

190 

0.046 

191 

0.326 

192 

0.436 

193 

0.026 

194 

65.117 

195 

65.118 

196 

34 .388 

197 

34.372 

198 

35.368 

199 

34.410 

200 

34.414 

201 

65.179 

202 

65.189 

203 

65 . 142 

204 

0.116 

205 

0.669 

206 

0.908 

207 

0.077 

208 

65.181 

209 

65.182 

210 

65.187 

211 

32.400 

212 

33.285 

213 

31.956 

214 

31.881 

215 

65 . 195 

216 

65.197 

217 

0.146 

218 

0.723 

219 

1.105 

220 

1.738 

221 

0.088 

222 

65.196 

223 

65.196 

224 

65.197 

225 

65.199 

226 

25.195 

227 

25.330 

228 

25,488 

229 

65.199 

230 

0.171 

231 

0.103 

232 

0.748 

233 

1.343 

234 

1.880 

235 

2.607 

236 

65.199 

237 

65.199 

238 

65.199 

239 

65.200 

240 

65.200 

241 

18.284 

242 

18.273 

243 

18.284 

244 

0.198 

245 

0.120 

246 

0.668 

247 

1.226 

248 

1.761 

249 

2.101 

250 

4.229 

251 

65.200 

252 

65.200 

253 

65.200 

254 

65.200 

255 

65.200 

256 

11.381 

257 

10.888 

258 

10 . 453 

259 

0.210 

260 

0.493 

261 

0.132 

262 

0.839 

263 

1.381 

264 

2.072 

265 

3.677 

266 

65.200 

267 

65.200 

268 

65.200 

“269 

65.200 

270 

65.200 

271 

2.559 

272 

0.221 

273 

0.519 

274 

0.862 

275 

0.134 

276 

0.954 

277 

0.938 

278 

1.950 

279 

65.200 

280 

65.200 

2 B 1 

65.200 

282 

65.200 

283 

65.200 

284 

1.635 

285 

1.394 

286 

0.204 

287 

0.501 

288 

0.814 

289 

0.900 

290 

0.130 

291 

65.200 

292 

65.200 

293 

65.200 

294 

65.200 

295 

65.200 

296 

1.128 

297 

1.100 

298 

0.158 

299 

0.379 

300 

0.670 

301 

0.605 

302 

0.099 

303 

65.200 

304 

65.200 

305 

65.200 

306 

65.200 

307 

65.200 

308 

0.913 

309 

0.094 

310 

0.060 

311 

0.240 

312 

0.444 

313 

0.397 

314 

65.200 

315 

65.200 

316 

65.200 

317 

65.200 

318 

65.200 

319 

0.044 

320 

0.028 

321 

0.103 

322 

0.143 

323 

0.145 

324 

65.200 

325 

65.200 

326 

65.200 

327 

65.200 

328 

65.200 

329 

0.011 

330 

0.028 

331 

0.006 

332 

0.036 

333 

0.035 

334 

65.200 

335 

65.200 

336 

65.200 

337 

65.200 

338 

65.200 

339 

0.002 

340 

0.005 

341 

0.001 

342 

0.006 

343 

0.006 

344 

65.200 

345 

65.200 

346 

65 . 200 

347 

65.200 

348 

65.200 

349 

0.000 

350 

0.001 

351 

0,001 

352 

0.000 

353 

0,001 

354 

65.200 

355 

65.200 

356 

65.200 

357 

65.200 

358 

65.200 

359 

0.000 

360 

0.000 

361 

0.000 

362 

0.000 

363 

- 0.000 

364 

65.200 

365 

65.200 

366 

65.200 

367 

65.200 

368 

65 . 200 

369 

0.000 

370 

0.000 

371 

0.000 

372 

0.000 

373 

0.000 

374 

65.200 

375 

65.200 

376 

65.200 

377 

65.200 

378 

65,200 

379 

- 0.000 

380 

0.000 

381 

0.000 

382 

0.000 

383 

- 0.000 

384 

65.200 

385 

65.200 

386 

65.200 

387 

65.200 

388 

65.200 

389 

0.000 

390 

0.000 

391 

0.000 

392 

0.000 

393 

0.000 

394 

65,200 

395 

65.200 

396 

65.200 

397 

65.200 

398 

65.200 

399 

- 0.000 

400 

0.000 

401 

0.000 

402 

0.000 

403 

- 0.000 

404 

65.200 
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Vl 

NEL 

S-ZETA 

S-EATA 

TUE-ZET-ET 


1 

72,753 

- 25.031 

0.097 

2 

128.730 

21.610 

- 10.909 

3 

150.471 

17.487 

- 3,885 

4 

134.717 

7.048 

- 3.030 

5 

103.118 

- 6.019 

- 6.335 

6 

42,985 

1.648 

1.088 

7 

43.985 

8.558 

0.592 

8 

40.778 

4.444 

0.450 

9 

36.630 

- 2.692 

0.738 

10 

34. 134 

- 8.519 

- 4.441 

11 

30.389 

- 9.585 

6.986 

12 

58.134 

- 24.751 

- 12.546 

13 

71.830 

- 7.608 

4.800 

14 

84.303 

- 0.846 

1.810 

15 

85.364 

4.561 

- 1.291 

16 

74.380 

5.382 

- 0.979 

17 

57.161 

0.881 

“ 0.938 

18 

182.210 

6.823 

- 4.466 

19 

193.090 

13.623 

- 6.950 

20 

184.950 

17.615 

- 4,893 

21 

138.810 

5.602 

2.921 

22 

141.180 

- 6.364 

3.446 

23 

159.040 

- 0.724 

8.298 

24 

159.540 

- 1.576 

- 4.006 

25 

143.690 

- 5 . 850 

- 3.097 

26 

181.040 

- 1.726 

- 2.301 

27 

175.540 

1.292 

- 2.021 

28 

206.120 

- 2.201 

- 0.493 

29 

196.080 

- 0.605 

- 2.915 

30 

219.495 

2.378 

8.540 

31 

220.655 

3.084 

3.971 

32 

219.076 

0.208 

- 4.697 

33 

218.126 

- 1.217 

- 8.853 

34 

220.024 

- 0.417 

2.685 

35 

219-586 

1.129 

3.916 

36 

219.361 

0.641 

0.488 

37 

216.837 

- 0.101 

- 3.518 

38 

212.754 

62.153 

24.500 

39 

281.875 

95.675 

- 25.342 

40 

22 S .291 

42.228 

- 19.808 

41 

198.368 

4.434 

- 18.915 

42 

196.845 

- 7.460 

11.578 

43 

205.437 

1.011 

6 . 168 

44 

211.376 

0.407 

- 3.798 

45 

218.019 

- 1.705 

- 8.118 

46 

54.648 

- 2.145 

- 4.063 

47 

60.289 

- 0.573 

“ 6.099 
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OF POOP. Q 



48 


64.369 

0.498 

-9.703 

49 


70.544 

2.808 

-12.586 

50 


71.693 

-12.204 

-2.760 

51 


44.375 

-42.569 

-23.700 

52 


-14.116 

-79.770 

-49.946 

53 


21.449 

105.675 

-27.785 

54 


“0.499 

-1.587 

-1.159 

55 


-9.181 

-0.671 

-2.958 

56 


18.219 

0.743 

-0.578 

57 


25.410 

0.581 

0.050 

58 


34,978 

1.825 

-0.734 

59 


40.194 

1.519 

-1.884 

60 


39.774 

1.679 

-5.518 

61 


32.715 

-3.919 

-8.640 

62 


-9.067 

-1.591 

1.703 

63 


-6.920 

“3.156 

1.557 

64 


-13.227 

-0.119 

0.270 

65 


-11.307 

-0.518 

1.700 

66 


-9.902 

2.126 

-1.099 

67 


-11.378 

0.579 

-2.753 

‘68 


-33.096 

-5.307 

3.648 

69 


-15.536 

1.261 

2.706 

70 


-36.337 

6.309 

-2.919 

71 


-40.329 

0.062 

-1.569 

72 


-44.515 

-2.928 

-0.573 

73 


-47.240 

-4.734 

-0.251 

74 


-55.141 

0.136 

-3.032 

75 


-55,361 

-1.528 

-1.849 

76 


-56.541 

-0.846 

“0,927 

77 


-58.733 

0.498 

-0.432 

78 


-62.359 

1.659 

0.193 

79 


-67.046 

-0.531 

5.309 

80 


-32.720 

-6.246 

-3.719 

81 


-38.555 

-5.031 

0.494 

82 


-44.718 

-3.183 

0.305 

83 


-50.290 

-0.961 

0.029 

84 


-54.344 

2.025 

0.268 

85 


-56.160 

4.501 

1.446 

86 


1.913 

2.701 

-3.104 

87 


1.889 

1.905 

-0.896 

88 


1.010 

1.776 

0.904 

89 


-5.577 

0.805 

3.054 

90 


1.411 

-4.648 

-8.905 



